% 
% hw5.m computes the 
% answers for hw5
%

clear all
diary hw5.mout
diary off
delete hw5.mout
diary hw5.mout

% ------ compute statistics for data using formulas ------ %
% - 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)' ] ;

% - enter formulas - %

display('answers from data')

g0data = dlyap1(A,Sigma);
g1data = A*g0data;
g2data = A^2*g0data;
g3data = A^3*g0data;

stdc = sqrt( g0data(2,2) )
stdn = sqrt( g0data(3,3) )

corc1c0 = g1data(2,2)/(stdc*stdc)
corc2c0 = g2data(2,2)/(stdc*stdc)
corc3c0 = g3data(2,2)/(stdc*stdc)

corrn0c = g0data(3,2)/(stdn*stdc)

% ------ compute statistics for model using formulas ------ %

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 ;
    0.2459  0.8858 ] ;

C = [0.2;0] ;

% - enter formulas - %

display('answers from model')

g0 = dlyap1(f2,C*C');
g1 = f2*g0;
g2 = f2^2*g0;
g3 = f2^3*g0;

varX0   = f1*g0*f1';
varX1X0 = f1*g1*f1';
varX2X0 = f1*g2*f1';
varX3X0 = f1*g3*f1';

stdc = sqrt( varX0(1,1) )
stdn = sqrt( varX0(2,2) )

corc1c0 = varX1X0(1,1)/(stdc*stdc)
corc2c0 = varX2X0(1,1)/(stdc*stdc)
corc3c0 = varX3X0(1,1)/(stdc*stdc)

corrn0c = varX0(2,1)/(stdn*stdc)

diary off