clear all;

% number of simulations
N=200;

Theta_s=zeros(N,13);
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.2;
variance_normal(4)=0.1;	
variance_normal(5)=0.1;
variance_normal(6)=0.05;
variance_normal(7)=0.09;
variance_normal(8)=0.05;
variance_normal(9)=0.05;
variance_normal(10)=0.009;
variance_normal(11)=0.009;
variance_normal(12)=0.009;
variance_normal(13)=0.009;

variance_normal=variance_normal/30;

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


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

% sorting the data in the right way
dy(:,1)=rw;
dy(:,2)=i;
dy(:,3)=dp;
dy(:,4)=y;

% Initial Values for parameters
load initialvalueGGLS.mat

% evaluating the first likelihood and prior functions
[loglike]=likeliggls(dy,Theta);
logpri=priorggls(Theta);

loglikelogpri=loglike+logpri;

save initiallikepriGGLS loglikelogpri

pau=0;

tic;
for j=1:N,
    
        newTheta=Theta+normrnd(0,variance_normal,1,13);
        newlogpri=priorggls(newTheta);
        if newlogpri>-Inf
                [newloglike]=likeliggls(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;

end

pau=pau/N*100;

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