% 
% ch7ex1.m simulates the economy from
% mid term 1 then estimates the 
% parameters of the model by GMM
%

clear all
diary ch7ex1.mout
diary off
delete ch7ex1.mout
diary ch7ex1.mout

% ----- set random number generator seed ----- %

randn('state',0);

% ------ set simulation length ----- %

t = 2000;

t2 = 2000;
t1 = 1999;
t0 = 1998;

d2 = 1000;
d1 = 999;
d0 = 998;

% ----- actual parameters used to solve economy ----- %

delta    = 0.025;
alpha    = 0.30;
thetabar = 1.05;
beta     = 0.99;
psi      = 1;
sigma    = 2;
phi      = 0.90;
sigmae   = 0.20;

% ----- solutions from midterm ----- %

f1 =[0.4681    0.4957
    0.3569   -0.1313]

f2 =[0.9000         0
    0.2459    0.8858 ]

C = [0.2;0]

% ----- compute states ----- % 

x(1:2,1)=[0;0];

for i = 1:t
 x(:,i+1) = f2*x(:,i) + C*randn(1);
end

% - compute endo vars - %

for i = 1:t
 y(:,i) = f1*x(:,i) ;
end

% ----- other endo vars go here ----- %
 
for i = 1:t
 out(:,i)    = .3*x(2,i) + (1-.3)*y(2,i) + x(1,i);
 invest(:,i) = (exp(x(2,i+1)) - (1-delta)*exp(x(2,i))/thetabar) ;
end

% ----- steady states from ch2ex1 ----- %
 
sol1 = [2.3501 0.4708 0.3731];
kbar = [2.3501];
cbar = [0.4708];
nbar = [0.3731];
ybar = kbar^.30*nbar^(1-.30) ;

% ----- compute ratios and re-define variables ----- %

ky2 = exp(x(2,d2:t2))*kbar./(exp(out(1,d2:t2))*ybar) ;
ky1 = exp(x(2,d1:t1))*kbar./(exp(out(1,d1:t1))*ybar) ;

cy2 = exp(y(1,d2:t2))*cbar./(exp(out(1,d2:t2))*ybar) ;
cy1 = exp(y(1,d1:t1))*cbar./(exp(out(1,d1:t1))*ybar) ;

gc2 = (exp(y(1,d2:t2))*cbar)./(exp(y(1,d1:t1))*cbar) ;
gc1 = (exp(y(1,d1:t1))*cbar)./(exp(y(1,d0:t0))*cbar) ;

n2  = exp(y(2,d2:t2))*nbar ;
n1  = exp(y(2,d1:t1))*nbar ;

lnz2 = x(1,d2:t2) ;
lnz1 = x(1,d1:t1) ;

% ----- now use csminwel.m to find the minimum 
%       of the squared moment conditions ----- %

omega = eye(4);

[fhat,xhat,ghat,hhat]=csminwel('ch7ex1F',[.80;.10;.40;1],eye(4),[],0.00001,1000,...
                                 lnz1,lnz2,gc1,gc2,ky1,ky2,n1,n2,cy1,cy2,omega,thetabar,beta,delta,sigma);
xhat

% ----- compute new optimal weighing matrix ----- %

ss = ch7ex2F(xhat,lnz1,lnz2,gc1,gc2,ky1,ky2,n1,n2,cy1,cy2,omega,thetabar,beta,delta,sigma);
omega = pinv(ss);

% ----- reoptimze gmm criterion with new omega ----- %

pause(2)
[fhat,xhat,ghat,hhat]=csminwel('ch7ex1F',xhat,eye(4),[],0.00001,1000,...
                                 lnz1,lnz2,gc1,gc2,ky1,ky2,n1,n2,cy1,cy2,omega,thetabar,beta,delta,sigma);
xhat

% ----- compute new optimal weighing matrix ----- %

ss = ch7ex2F(xhat,lnz1,lnz2,gc1,gc2,ky1,ky2,n1,n2,cy1,cy2,omega,thetabar,beta,delta,sigma);
omega = pinv(ss)