%
% hw3sol.m is an example matlab program 
%

% ----- save output to file ----- %

clear all
diary hw3sol.mout
diary off
delete hw3sol.mout
diary hw3sol.mout

% ----- display date and time of computation ----- %

format bank
date
time0 = clock;

% ----- compute first and second moments ----- %
% ------ ex2.4(bi) ----- %

rho1  = 1.2;
rho2  = -0.3;
rho3  = 0;
rho4  = 0;
mu    = 10;
c     = 1;
alpha = mu*(1-rho1 - rho2 - rho3 - rho4);

A = [ rho1 rho2 rho3 rho4 alpha; 
      1 0 0 0 0;
      0 1 0 0 0;
      0 0 1 0 0;
      0 0 0 0 1 ] ;

C = [c; 0; 0; 0; 0] ;

% - covariance stat mean 
%   is associated with the 
%   unit eigenvalue in D -%

[mux D] = eig(A) 
Cx      = doublej(A,C*C')

% ------ ex2.4(bii) ------ %

rho1  = 1.2;
rho2  = -0.3;
rho3  = 0;
rho4  = 0;
c     = 2;
mu    = 10;
alpha = mu*(1-rho1 - rho2 - rho3 - rho4);

A = [ rho1 rho2 rho3 rho4 alpha; 
      1 0 0 0 0;
      0 1 0 0 0;
      0 0 1 0 0;
      0 0 0 0 1 ] ;

C = [c; 0; 0; 0; 0] ;

% - covariance stat mean 
%   is associated with the 
%   unit eigenvalue in D -%

[mux D] = eig(A) 
Cx      = doublej(A,C*C')


% ------ ex2.4(biii) ----- %

rho1  = 0.90;
rho2  = 0;
rho3  = 0;
rho4  = 0;
mu    = 5;
c     = 1;
alpha = mu*(1-rho1 - rho2 - rho3 - rho4);

A = [ rho1 rho2 rho3 rho4 alpha; 
      1 0 0 0 0;
      0 1 0 0 0;
      0 0 1 0 0;
      0 0 0 0 1 ] ;

C = [c; 0; 0; 0; 0] ;

% - covariance stat mean 
%   is associated with the 
%   unit eigenvalue in D -%

[mux D] = eig(A) 
Cx      = doublej(A,C*C')

% ------ ex2.4(biv) ----- %

% ------ ex2.4(bv) ----- %

% ------ ex2.4(c) ----- %
% - we know that Et[y(t+5)] is the first row and column of A0^5*x(t) 
%   thus we need to compute A0^5 and multiply x(t) by it to find the h's - %

% ------ ex2.4(ci) ----- %

rho1  = 1.2;
rho2  = -0.3;
rho3  = 0;
rho4  = 0;
mu    = 10;
c     = 1;
alpha = mu*(1-rho1 - rho2 - rho3 - rho4);

A = [ rho1 rho2 rho3 rho4 alpha; 
      1 0 0 0 0;
      0 1 0 0 0;
      0 0 1 0 0;
      0 0 0 0 1 ] ;

C = [c; 0; 0; 0; 0] ;

A^5

% - thus we see 0.7387*y(t)-0.2603*y(t-1)+0*y(t-2)+0*y(t-3)+5.2156*1  - %

% ------ ex2.4(cii) ----- %

% ------ ex2.4(ciii) ----- %

% ------ ex2.4(civ) ----- %

% ------ ex2.4(cv) ----- %

% ------ ex2.4(d) ----- %
% - we need to compute  Et[sum{.95^j*y(t+j)}]. Note that this is essentially
%   Et[sum{.95^j*[1 0 0 0 0]*x(t+j)}]. We know how to compute this conditional
%   expectation; sum{.95^j*[1 0 0 0 0]*A^jx(t)} = [1 0 0 0 0]*inv(I -beta*A)*x(t).  - %

% ------ ex2.4(di) ----- %

rho1  = 1.2;
rho2  = -0.3;
rho3  = 0;
rho4  = 0;
mu    = 10;
c     = 1;
alpha = mu*(1-rho1 - rho2 - rho3 - rho4);

A = [ rho1 rho2 rho3 rho4 alpha; 
      1 0 0 0 0;
      0 1 0 0 0;
      0 0 1 0 0;
      0 0 0 0 1 ] ;

C = [c; 0; 0; 0; 0] ;

[1 0 0 0 0]*inv(eye(5) - 0.95*A)

% - thus we see 7.6482*y(t)-2.179*y(t-1)+0*y(t-2)+0*y(t-3)+145.3155*1  - %

% ------ ex2.4(dii) ----- %

% ------ ex2.4(diii) ----- %

% ------ ex2.4(div) ----- %

% ------ ex2.4(dv) ----- %

% ------ ex2.4(e) ----- %
% - we need to compute  Cx(j) = A^j*Cx(0) - %

% ------ ex2.4(ei) ----- %

rho1  = 1.2;
rho2  = -0.3;
rho3  = 0;
rho4  = 0;
mu    = 10;
c     = 1;
alpha = mu*(1-rho1 - rho2 - rho3 - rho4);

A = [ rho1 rho2 rho3 rho4 alpha; 
      1 0 0 0 0;
      0 1 0 0 0;
      0 0 1 0 0;
      0 0 0 0 1 ] ;

C = [c; 0; 0; 0; 0] ;


Cx0  = doublej(A,C*C')

Cx1  = A*Cx0 ;
Cx5  = A^5*Cx0 ;
Cx10 = A^10*Cx0 ;


Cx1(1,1)
Cx5(1,1)
Cx10(1,1)


% - we see that Cx1=6.85, Cx5=3.70, Cx10=1.59  - %

% ------ ex2.4(eii) ----- %

% ------ ex2.4(eiii) ----- %

% ------ ex2.4(eiv) ----- %

% ------ ex2.4(ev) ----- %


% ----- turn off diary ----- %

comptime = etime(clock, time0)
diary off