% 
% ch4ex1.m simulates the economy from
% mid term 1
%


clear all
diary ch4ex1.mout
diary off
delete ch4ex1.mout
diary ch4ex1.mout




% ----- simulate the economy from midterm ----- %

delta = 0.025;
alpha = 0.30;
thetabar = 1.05;

%- 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:1000
 x(:,i+1) = f2*x(:,i) + C*randn(1);
end

% - compute endo vars - %

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

% - other endo vars go here - %
 
for i = 1:1000
 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

% ----- compute sample statistics for economy ---- %

std([out(1,100:1000)' y(1,100:1000)' y(2,100:1000)' invest(1,100:1000)'])
corrcoef([out(1,100:1000)' y(1,100:1000)' y(2,100:1000)' invest(1,100:1000)'])

% ---- VAR representation of sytem ----- %

A1=[1 0 0 0; 
    0 1 0 0; 
    -0.4681 -0.4957 1 0; 
    -0.3569 0.1313 0 1];

A2 = [0.90 0 0 0; 
      0.2459 .8858 0 0; 
      0 0 0 0; 
      0 0 0 0];
  
pinv(A1)*A2
pinv(A1)

% ---- alternative VAR representation of sytem ----- %

A1=[1 -0.30 0 -0.70; 
    0 1 0 0; 
    -0.4681 (0.4681*.30-0.4957) 1 (0.4681*0.70); 
    -0.3569 (0.3569*.30+0.1313) 0 (1+0.3569*0.70)];

A2 = [0.90 (-0.90*0.30) 0 (-0.90*0.70); 
      0.2459 (-0.2459*0.30+.8858) 0 (-0.2459*.70); 
      0 0 0 0; 
      0 0 0 0];
  
pinv(A1)*A2
pinv(A1)

% ----- now using fake data, run var to 
%       get covariance structure of economy ----- %

% - drop first 100 and redefine data and lags - %

top = 1000
top1 = 999
drop = 800
drop1 = 799

gdp  = out(1,drop:top)';
con  = y(1,drop:top)';
hour = y(2,drop:top)'; 
inve = invest(1,drop:top)';

L1gdp  = out(1,drop1:top1)';
L1con  = y(1,drop1:top1)';
L1hour = y(2,drop1:top1)' ;
L1inve = invest(1,drop1:top1)';

% ----- now use regression formula ----- %

Y1 = [gdp];
Y2 = [con];
Y3 = [hour];
Y4 = [inve];

Y = [Y1;Y2;Y3;Y4] ;

X1 = [ones(top-drop+1,1) L1gdp L1con L1hour L1inve];
X = kron(eye(4),X1);

betahat = pinv(X'*X)*X'*Y 

% ---- compute residuals and thus var/covar matrix of innovations - %

u1 = Y1 - X1*betahat(1:5);
u2 = Y2 - X1*betahat(6:10);
u3 = Y3 - X1*betahat(11:15);
u4 = Y4 - X1*betahat(16:20);

Sigma = 1/(top-drop+1)*[u1'*u1 u1'*u2 u1'*u3 u1'*u4 ;
                        u2'*u1 u2'*u2 u2'*u3 u2'*u4 ;
                        u3'*u1 u3'*u2 u3'*u3 u3'*u4 ;
                        u4'*u1 u4'*u2 u4'*u3 u4'*u4 ]
           
% ----- compute covariances ----- %

% - reshape beta - %

A = [ betahat(2:5)',
      betahat(7:10)' ,
      betahat(12:15)' ,
      betahat(17:20)' ] ;
           
Gamma0 = dlyap1(A,Sigma) 
Gamma1 = A*Gamma0 
Gamma2 = A^2*Gamma0 
Gamma3 = A^3*Gamma0 

% ----- compute correlations ----- %

corr=[Gamma0(1,2)/(sqrt(Gamma0(1,1))*sqrt(Gamma0(2,2)));
     Gamma0(1,3)/(sqrt(Gamma0(1,1))*sqrt(Gamma0(3,3)));
     Gamma0(1,4)/(sqrt(Gamma0(1,1))*sqrt(Gamma0(4,4)))]

% ----- compute autocorrelations ----- %

autocorryc=[Gamma0(1,2)/(sqrt(Gamma0(1,1))*sqrt(Gamma0(2,2)));
            Gamma1(1,2)/(sqrt(Gamma0(1,1))*sqrt(Gamma0(2,2)));
            Gamma2(1,2)/(sqrt(Gamma0(1,1))*sqrt(Gamma0(2,2)));
            Gamma3(1,2)/(sqrt(Gamma0(1,1))*sqrt(Gamma0(2,2)))]

autocorryi=[Gamma0(1,4)/(sqrt(Gamma0(1,1))*sqrt(Gamma0(4,4)));
            Gamma1(1,4)/(sqrt(Gamma0(1,1))*sqrt(Gamma0(4,4)));
            Gamma2(1,4)/(sqrt(Gamma0(1,1))*sqrt(Gamma0(4,4)));
            Gamma3(1,4)/(sqrt(Gamma0(1,1))*sqrt(Gamma0(4,4)))]

% ----- data's impulse responses using no theory  ----- %
        
% - start to compute recursive impulse - %

 An(:,:,1) = eye(4);

  iters = 50;
 
 for i = 1:iters
  An(:,:,i+1) = A*An(:,:,i);  
 end
 clear i
 
 for i = 1:4
  for j = 1:iters
  chon(:,j,i) = Sigma(i,i)^(-1/2)*An(:,:,j)*Sigma*An(:,i,1);
  end 
  clear j
 end
 clear i
 
% ----- compute models impulses ----- %

 chonm = C;
 
 for i = 1:iters
 chonm = [chonm f1*f2^i*C ];
end
 
figure(1)
 subplot(2,1,1)
 plot([chon(2:3,:,1)'])
 title('data impulses for c and n')
 axis([0 iters -0.05 .15])
 
 subplot(2,1,2)
 plot([chonm'])
 title('model impulses for c and n')
 axis([0 iters -0.05 .15])

% ----- read in actual data ----- %

X = load('rbc1.txt') ;

gdp  = log(X(:,1));
con  = log(X(:,2));
hour = log(X(:,4));
inve = log(X(:,3));

% ----- hp filter data ----- %


tgdp  = hpfilter(gdp,1600);
tcon  = hpfilter(con,1600);
thour = hpfilter(hour,1600);
tinve = hpfilter(inve,1600);

sgdp  = gdp - tgdp;
scon  = con-tcon;
shour = hour-thour;
sinve = inve-tinve;

% ---- run regression ------ %

[s1,s2] = size(sgdp);

ss1 = s1-1;

Y1 = [sgdp(2:s1,1)];
Y2 = [scon(2:s1,1)];
Y3 = [shour(2:s1,1)];
Y4 = [sinve(2:s1,1)];

Y = [Y1;Y2;Y3;Y4] ;

X1 = [ones(s1-1,1) sgdp(1:ss1,1) scon(1:ss1,1) shour(1:ss1,1) sinve(1:ss1,1)];
X = kron(eye(4),X1);

betahat = pinv(X'*X)*X'*Y 

% ---- estimate sigma ------ %

u1 = Y1 - X1*betahat(1:5);
u2 = Y2 - X1*betahat(6:10);
u3 = Y3 - X1*betahat(11:15);
u4 = Y4 - X1*betahat(16:20);

Sigma = 1/(ss1-5)*[u1'*u1 u1'*u2 u1'*u3 u1'*u4 ;
                        u2'*u1 u2'*u2 u2'*u3 u2'*u4 ;
                        u3'*u1 u3'*u2 u3'*u3 u3'*u4 ;
                        u4'*u1 u4'*u2 u4'*u3 u4'*u4 ]
           

% - reshape beta - %

A = [ betahat(2:5)',
      betahat(7:10)' ,
      betahat(12:15)' ,
      betahat(17:20)' ] ;

Gamma0 = dlyap1(A,Sigma) 
Gamma1 = A*Gamma0 
Gamma2 = A^2*Gamma0 
Gamma3 = A^3*Gamma0