clear all;

% adjustment parameter

% parameters of the simulation

%were are saving the # of simulations 

N=2000; % number of simulations
k=8; %number of lags

Theta_s=zeros(N,15);
loglike_s=zeros(N,1);
logpri_s=zeros(N,1);

% variance of metropolis-hasting 

variance_normal(1)=0.16;
variance_normal(2)=0.1;
variance_normal(3)=0.1;
variance_normal(4)=0.2;
variance_normal(5)=0.2;
variance_normal(6)=0;
variance_normal(7)=0.1;
variance_normal(8)=0.05;
variance_normal(9)=0.09;
variance_normal(10)=0.05;
variance_normal(11)=0.05;
variance_normal(12)=0.009;
variance_normal(13)=0.009;
variance_normal(14)=0.009;
variance_normal(15)=0.009;

variance_normal=variance_normal/20;

% reading the data

[rw,i,dp,y]=textread('usadefl1d.txt','%f %f %f %f');

dy=zeros(size(dp,1),4); %data

dy(:,1)=rw;
dy(:,2)=i;
dy(:,3)=dp;
dy(:,4)=y;

% Initial Values for parameters

load initialvalueEHL.mat

% evaluating the first likelihood and prior functions

[loglike]=likeliehl(dy,Theta);
logpri=priorehl(Theta);

loglikelogpri=loglike+logpri;

save initiallikepriEHL loglikelogpri

% initial value for the acceptance rate, pau is acceptance rate in honor to my
% lovely co-author

pau=0;
test1=0;

%fid10=fopen('thetaehl.txt','w');
%fid11=fopen('likepri.txt','w');
%fid12=fopen('ACPI.txt','w');
%fid13=fopen('ACY.txt','w');
%fid14=fopen('ACRW.txt','w');
%fid15=fopen('ACI.txt','w');
%fid16=fopen('CCYPI.txt','w');
%fid17=fopen('CCYRW.txt','w');

tic;
for j=1:N,

        newTheta=Theta+normrnd(0,variance_normal,1,15);
	    newTheta(5)=4.33;
        newTheta(6)=0;                        %Set the rule of thumb behavior to zero^M
        newlogpri=priorehl(newTheta);
        if newlogpri>-Inf
                [newloglike]=likeliehl(dy,newTheta);
                if newloglike==-Inf;
                        ratio=0;
                else
                        ratio=exp(newloglike+newlogpri-(loglike+logpri));
                        if rand<=ratio % we accept with prob ratio
                                loglike=newloglike;
                                logpri=newlogpri;
                                Theta=newTheta;
                                pau=pau+1;
                        end
                end
        end

	Theta_s(j,:)=Theta;
	loglike_s(j,1)=loglike;
	logpri_s(j,1)=logpri;

   %	ACPI=crosscorrpi(Theta,k);
   %	[ACY,CCYPI,CCYRW]=crosscorry(Theta,k);
   %	ACRW=crosscorrrw(Theta,k);
   %	ACI=crosscorri(Theta,k);
    
end
toc;

% saving rate of acceptance

pau=pau/N*100;
test1=test1/N*100;

fid4=fopen('EHLrate.txt','w');
fprintf(fid4,'%s ','# of draws:');
fprintf(fid4,'%d\n',N);
fprintf(fid4,'%s ','rate of acceptance (in %):');
fprintf(fid4,'%f\n',pau);
fprintf(fid4,'%s ','% of indeterminancies:');
fprintf(fid4,'%f\n',test1);   
fclose(fid4);

save EHLoutput