% P E M _ e s t i m _ U S
% ***********************

 close all
 clear all

format 
set(0,'DefaultLineLineWidth',1);
set(0,'DefaultAxesFontSize',8);

% Loading of data
 load USAMACRO;

% Variables
% #########

  usamacro = USAMACRO;
  
 [N,n2]=size(usamacro(:,:));
% t
 C   = usamacro(4:N-1,2);  % C(t)
 I   = usamacro(4:N-1,3);  % I(t)
 R   = usamacro(4:N-1,4);  % R(t)
 Y   = usamacro(4:N-1,5);  % Y(t)
 G   = usamacro(4:N-1,6);  % G(t)
 M   = usamacro(4:N-1,7);  % M(t)
% t-1
 Y_1 = usamacro(3:N-2,5);  % Y(t-1)
 R_1 = usamacro(3:N-2,4);  % R(t-1)
% t-2
 R_2 = usamacro(2:N-3,4);  % R(t-2)
% t-3
 R_3 = usamacro(1:N-4,4);  % R(t-3)
% dt
 dY  = Y - Y_1;             % dY=Y(t)-Y(t-1)
% t+1
 Cnew = usamacro(5:N,2);    % C(t+1)
 Inew = usamacro(5:N,3);    % I(t+1)
 Rnew = usamacro(5:N,4);    % R(t+1)
 Ynew = usamacro(5:N,5);    % Y(t+1)
 Gnew = usamacro(5:N,6);    % G(t+1)
 Mnew = usamacro(5:N,7);    % M(t+1)
% Difference
 dMnew = Mnew -M;            % dM(t+1)=M(t+1)-M(t)
 dYnew = Ynew - Y;           % dY(t+1)=Y(t+1)-Y(t)
% Constant
 const=ones(size(Cnew));  % constant
% 
 t=usamacro(5:N,1);

%...diary estim_us.txt %<---------<<<
disp('  ');
disp('  ');
disp('  ');
disp('   SMALL LINEAR MACROECONOMIC MODEL OF   U S A   ECONOMY');
disp('   #####################################################');
disp('   ');
disp('   ');
disp('   V a r i a b l e s');
disp('   C    real aggregate personal consumption    ');
disp('   I    real gross domestic investment         ');
disp('   Y    real GNP (net of exports and imports)  ');
disp('   G    real government spending               ');
disp('   M    real money stock (M1)                  ');
disp('   R    interest rate on 3-month Treasury bills');
disp('  ');
disp('  ');
disp('   R e g r e s s   E s t i m a t i o n  of Individual Equations');
disp('   ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++');
disp('  ');
disp('   Hypothese No.1: C_t+1= c0 * const  +  c1 * Y_t+1  +   c2 * C_t');
disp('   **************************************************************');

% Regress Estimation No.1
% ***********************
% Linearni regrese
% ****************  
% Vysvetlovana promenna
  y1=Cnew;
% Rozsirena matice soustavy ( matice planu )
  X=[const Ynew C];
% Pocet pozorovani a pocet parametru
  [n k]=size(X);
% Linearni regrese
  [b,b_int,e,e_int,STATS] = regress(y1,X,0.05);
  c = b;
  ec= e;
% Vypoctene hodnoty endogenni promenne
  y1r=X*b;

% Standardni chyby
% ****************
% Rozptyl chyby vyrovnani
  sigma2=(e'*e)/(n-k);
% Standardni odchylky parametru
  Sb=sigma2*(X'*X)^(-1);
  stderrb=sqrt(diag(Sb));

% t testy
% *******
% t testy parametru
  ttest = abs(b)./stderrb;
% kriticka hodnota t testu
  t_kr  = tinv(0.975,n-k);

%...............................................

disp('   Regresni koeficienty')   
disp('   co        c1        c2');
disp(b')
disp('   Standardni odchylky koeficientu');
disp('   std_c0    std_c1    std_c2');
disp(stderrb')
disp('   t-testy')
disp('   t_c0      t_c1      t_c2');
disp(ttest')
disp('   Kriticka hodnota t*');
disp(t_kr);
%...............................................

% Korigovany koeficient determinace
% *********************************
% Koeficient determinace
  R2 = STATS(1,1);
% Korigovany koeficient determinace
  R2A = 1-(1-R2)*((n-1)/(n-k));

% Durbin-Watsonova statistika
% ***************************
  dw1 = 0;
 for i=2:n
	  dw=(e(i)-e(i-1))^2;
        dw1=dw1+dw;
 end
  dw = dw1/sum(e.^2);
 dwc=dw;
%...........................................
disp('     RA^2       DW   krit.hodnota t');
disp([R2A  dw tinv(0.975,n-k) ]);
%...........................................

figure(1)
clf
set(gcf,'Name','Observed Values and Predicted Values of C  (Method: " O L S ")');
subplot(2,1,1)
plot(t,Cnew,'r');
hold on
plot(t,y1r,'b-.');
hold off
axis([min(t)-1,max(t)+1,min(Cnew)-60,max(Cnew)+60]);
title('Real Aggrecate Persona Consumption C');
legend('observed C','predicted C ( " O L S " )');
ylabel('C(t), C_p(t)   [ billions of 1982 dollars ]');
xlabel('Time   [ qurter ]');
grid
subplot(2,1,2)
plot(t,e,'r');
hold on
plot(t,[e_int(:,1)-e,e_int(:,2)-e],'b:');
hold off
axis([min(t)-1,max(t)+1,min(e)-5,max(e)+5]);
legend('errors','confidence interval');
title('Errorbars of Consumption  C');
ylabel('Errors of C(t)   [ billions of 1982 dollars ]');
xlabel('Time   [ qurter ]');
grid

disp('   Hypothese No.2: I_t+1 = i0  +  i1*dY_t  +  i2*Y_t+1 + i3*R_t-3');
disp('   **************************************************************');

% Regress Estimation No.2
% ***********************
% Linearni regrese
% ****************  
% Vysvetlovana promenna
  y2=Inew;
% Rozsirena matice soustavy ( matice planu )
  X2=[const dY Ynew R_3];
% Pocet pozorovani a pocet parametru
  [n k]=size(X2);
% Linearni regrese
  [b,b_int,e,e_int,STATS] = regress(y2,X2,0.05);
  ii = b;
  ei = e;
% Vypoctene hodnoty endogenni promenne
  y2r=X2*b;

% Standardni chyby
% ****************
% Rozptyl chyby vyrovnani
  sigma2=(e'*e)/(n-k);
% Standardni odchylky parametru
  Sb=sigma2*(X2'*X2)^(-1);
  stderrb=sqrt(diag(Sb));

% t testy
% *******
% t testy parametru
  ttest = abs(b)./stderrb;
% kriticka hodnota t testu
  t_kr  = tinv(0.975,n-k);

%...............................................

disp('   Regresni koeficienty')   
disp('   i0        i1        i2        i3');
disp(b')
disp('   Standardni odchylky koeficientu');
disp('   std_i0    std_i1    std_i2    std_i3');
disp(stderrb')
disp('   t-testy')
disp('   t_i0      t_i1      t_i2      t_i3');
disp(ttest')
disp('   Kriticka hodnota t*');
disp(t_kr);
%...............................................

% Korigovany koeficient determinace
% *********************************
% Koeficient determinace
  R2 = STATS(1,1);
% Korigovany koeficient determinace
  R2A = 1-(1-R2)*((n-1)/(n-k));

% Durbin-Watsonova statistika
% ***************************
  dw1 = 0;
 for i=2:n
	  dw=(e(i)-e(i-1))^2;
        dw1=dw1+dw;
 end
  dw = dw1/sum(e.^2);
 dwi=dw;
%...........................................
disp('     RA^2       DW   krit.hodnota t');
disp([R2A  dw tinv(0.975,n-k) ]);
%...........................................


figure(2)
clf
set(gcf,'Name','Observed Values and Predicted Values of I  ( Method: " O L S ")');
subplot(2,1,1)
plot(t,Inew,'r');
hold on
plot(t,y2r,'b-.');
hold off
axis([min(t)-1,max(t)+1,min(Inew)-30,max(Inew)+30]);
title('Real Gross Domestic Investment I');
legend('observed I','predicted I ( " O L S " )');
ylabel('I(t),  I_p(t)   [ billions of 1982 dollars ]');
xlabel('Time   [ quarter ]');
grid

subplot(2,1,2)
plot(t,e,'r');
hold on
plot(t,[e_int(:,1)-e,e_int(:,2)-e],'b:');
hold off
axis([min(t)-1,max(t)+1,min(e)-5,max(e)+5]);
legend('errors','confidence interval');
title('Errorbars of  I');
ylabel('Errors of I(t)   [ billions of 1982 dollars ]');
xlabel('Time   [ quarter ]');
grid

disp('   Hypothese No.3:  R_t+1 =  r0  +  r1*Y_t+1  +  r2*dY_t+1  ');
disp('                             + r3*dM_t+1  + r4*(R_t + R_t-1)');
disp('   *********************************************************'); 
% Regress Estimation No.3
% ***********************
% Linearni regrese
% ****************  
% Vysvetlovana promenna
  y3=Rnew;
% Rozsirena matice soustavy ( matice planu )
  X3=[const Ynew dYnew dMnew (R + R_1)];
% Pocet pozorovani a pocet parametru
  [n k]=size(X3);
% Linearni regrese
  [b,b_int,er,e_int,STATS] = regress(y3,X3,0.05);
  r = b;
% Vypoctene hodnoty endogenni promenne
  y3r=X3*b;

% Standardni chyby
% ****************
% Rozptyl chyby vyrovnani
  sigma2=(er'*er)/(n-k);
% Standardni odchylky parametru
  Sb=sigma2*(X3'*X3)^(-1);
  stderrb=sqrt(diag(Sb));

% t testy
% *******
% t testy parametru
  ttest = abs(b)./stderrb;
% kriticka hodnota t testu
  t_kr  = tinv(0.975,n-k);

%...............................................

disp('   Regresni koeficienty')   
disp('   r0        r1        r2        r3        r4');
disp(b')
disp('   Standardni odchylky koeficientu');
disp('   std_r0    std_r1    std_r2    std_r3    std_r4');
disp(stderrb')
disp('   t-testy')
disp('   t_r0      t_r1      t_r2      t_r3      t_r4');
disp(ttest')
disp('   Kriticka hodnota t*');
disp(t_kr);
%...............................................

% Korigovany koeficient determinace
% *********************************
% Koeficient determinace
  R2 = STATS(1,1);
% Korigovany koeficient determinace
  R2A = 1-(1-R2)*((n-1)/(n-k));

% Durbin-Watsonova statistika
% ***************************
  dw1 = 0;
 for i=2:n
	  dw=(e(i)-e(i-1))^2;
        dw1=dw1+dw;
 end
  dw = dw1/sum(e.^2);
 dwr=dw;
%...........................................
disp('     RA^2       DW   krit.hodnota t');
disp([R2A  dw tinv(0.975,n-k) ]);
%...........................................

figure(3)
clf
set(gcf,'Name','Obseved Values and Predicted Values of R  (Method: " O L S ")');
subplot(2,1,1)
plot(t,Rnew,'r');
hold on
plot(t,y3r,'b-.');
hold off
axis([min(t)-1,max(t)+1,min(Rnew)-0.5,max(Rnew)+0.5]);
title('Interest Rate on 3 Month Treasury Bills  R');
legend('observd R','predicted R ( " O L S " )');
ylabel('R(t),   R_p(t)   [ percent per year ]');
xlabel('Time  [ quarter ]');
grid
subplot(2,1,2)
plot(t,er,'r');
hold on
plot(t,[e_int(:,1)-er,e_int(:,2)-er],'b:');
hold off
axis([min(t)-1,max(t)+1,min(er)-0.5,max(er)+0.5]);
legend('errors','confidence interval');
title('Errorbars of R');
ylabel('Errors of  R(t)   [ percent per year ]');
xlabel('Time  [ quarter ]');
grid

figure(4)
clf
set(gcf,'Name','Obseved Values and Predicted Values of Y  (Method: " O L S ")');
subplot(2,1,1)
plot(t,Ynew,'r');
hold on
plot(t,y1r+y2r+Gnew,'b-.');
hold off
axis([min(t)-1,max(t)+1,min(y1r+y2r+Gnew)-160,max(y1r+y2r+Gnew)+160]);
title('Gross National Product');
legend('observd Y','predicted Y ( " O L S " )');
ylabel('Y(t),   Y_p(t)   [ billions of 1982 dollars ]');
xlabel('Time  [ quarter ]');
grid

subplot(2,1,2)
plot(t,Ynew-y1r-y2r-Gnew,'r');
hold off
axis([min(t)-1,max(t)+1,min(Ynew-y1r-y2r-Gnew)-10,max(Ynew-y1r-y2r-Gnew)+10]);
legend('errors','confidence interval');
title('Errorbars of Y');
ylabel('Errors of  Y(t)   [ billions of 1982 dollars ]');
xlabel('Time  [ quarter ]');
grid

% S I M U L T A N E O U S  S I M U L A T I O N  O F  S T A T E   M O D E L
% with parameters estimated by OLS on individual equations  
% ########################################################################
% 
% State space model
% *****************
% state equation
%                   x_t+1 = A x_t  +  B u_t + w(t)
% output equation
%                   y_t   = C x_t  +  D u_t + e(t)
% initial condition
%                   x_to  = x0

%  Original econometric equations
%  C_t+1 = c0 * const  +  c1 * Y_t+1  +   c2 * C_t                                ( 1 )
%  I_t+1 = i0 * const  +  i1*dY_t  +  i2*Y_t+1 + i3*R_t-3                         ( 2 )
%  R_t+1 = r0 * const  +  r1*Y_t+1  +  r2*dY_t+1 + r3*dM_t+1  + r4*(R_t + R_t-1)  ( 3 )
%  Y_t+1 = C_t+1  + I_t+1  + G_t+1                                                ( 4 )

% Aprior Estimated Parameters of Individual Equations
% ---------------------------------------------------
 c0 = c(1,1);  c1 = c(2,1);  c2 = c(3,1);
 i0 = ii(1,1); i1 = ii(2,1); i2 = ii(3,1); i3 = ii(4,1);
 r0 = r(1,1);  r1 = r(2,1);  r2 = r(3,1);  r3 = r(4,1); r4 = r(5,1);
 
% M a t r i c e s   o f   o r i g i n a l   e c o n o m e t r i c  m o d e l 
%  x_t+1 = Ao x_t+1 + A1 x_t + B1 u_t+1
% **************************************************************************

% M a t r i x   Ao
% ----------------
% 1     2     3     4     5     6     7     8      
% C_t+1 I_t+1 R_t+1 Y_t+1 R_t   Y_t   R_t-1 R_t-2   
% __________________________________________________________
Ao=[
  0     0     0     c1    0     0     0     0     % 1  C_t+1
  0     0     0     i2    0     0     0     0     % 2  I_t+1 
  0     0     0   (r1+r2) 0     0     0     0     % 3  R_t+1
  1     1     0     0     0     0     0     0     % 4  Y_t+1 
  0     0     0     0     0     0     0     0     % 5  Y_t 
  0     0     0     0     0     0     0     0     % 6  R_t
  0     0     0     0     0     0     0     0     % 7  R_t-1
  0     0     0     0     0     0     0     0];   % 8  R_t-2

% M a t r i x   A1
% ----------------
% 1     2     3     4     5     6     7     8     
% C_t   I_t   R_t   Y_t   Y_t-1 R_t-1 R_t-2 R_t-3
% __________________________________________________________
A1 = [
  c2    0     0     0     0     0     0     0     % 1  C_t+1
  0     0     0     i1   -i1    0     0     i3    % 2  I_t+1 
  0     0     r4   -r2    0     r4    0     0     % 3  R_t+1
  0     0     0     0     0     0     0     0     % 4  Y_t+1 
  0     0     0     1     0     0     0     0     % 6  Y_t 
  0     0     1     0     0     0     0     0     % 6  R_t
  0     0     0     0     0     1     0     0     % 7  R_t-1
  0     0     0     0     0     0     1     0];   % 8  R_t-2

% M a t r i x   B1
% ----------------
% 1      2     3
% const, Gnew, dMnew 
% _____________________________
B1 = [
  c0     0     0     % 1  C_t+1
  i0     0     0     % 2  I_t+1 
  r0     0     r3    % 3  R_t+1
  0      1     0     % 4  Y_t+1 
  0      0     0     % 5  Y_t   
  0      0     0     % 6  R_t 
  0      0     0     % 7 R_t-1
  0      0     0];   % 8 R_t-2
% Matrix C
% --------
 CC = zeros(4,8);

%  1    2     3    4
%  C    I     R    Y
%  _________________
 CC(1:4,1:4)=[
   1    0     0    0  % 1 C
   0    1     0    0  % 2 I
   0    0     1    0  % 3 R  
   0    0     0    1];% 4 Y

% Matrix D
% --------
 D=zeros(4,3);

%  1     2   3        
%  const G   dMnew   
%_____________________________
 D(1:4,1:3)=[
   0     0   0   % 1 C
   0     0   0   % 2 I 
   0     0   0   % 3 R
   0     0   0]; % 4 Y

% Matrix K
% --------
 K=zeros(8,4);
 K(1,1)=0;
 K(2,2)=0;
 K(3,3)=0;
 K(4,4)=0;
 
% S t a t e   m o d e l
% #################################################  
% x(t+1)  = A x(t)  +  B u(t+1)  +  K e_t
%   y(t)  = C x(t)  +  e_t 
%   x(0)  = x0
% _________________________________________________
% A(t) = {[ I - Ao(t) ]^(-1)}*A1(t)
% B(t) = {[ I - Ao(t) ]^(-1)}*B1(t)
% C(t) = {[ I - Ao(t) ]^(-1)}*C1(t)
% Ao, A1, B1, C1 jsou matice ekonometrickeho modelu
% #################################################
  n = 8; % number of states
  m = 4; % number of outputs
  r = 3; % number of inputs

% Restructured Matrices of State Model
% ************************************
  Io = eye(n);

  A = inv([Io-Ao])*A1;
  
  B = inv([Io-Ao])*B1;
  
% New Matrices of State Model
% ***************************

% Matrix C
% --------
 CC = zeros(4,8);

%  1    2     3    4
%  C    I     R    Y
%  _________________
 CC(1:4,1:4)=[
   1    0     0    0  % 1 C
   0    1     0    0  % 2 I
   0    0     1    0  % 3 R  
   0    0     0    1];% 4 Y

% Matrix D
% --------
 D=zeros(4,3);

%  1     2   3        
%  const G   dMnew   
%_____________________________
 D(1:4,1:3)=[
   0     0   0   % 1 C
   0     0   0   % 2 I 
   0     0   0   % 3 R
   0     0   0]; % 4 Y

% Matrix K
% --------
 K=zeros(8,4);
 K(1:4,1:4)=eye(4);
 
% Outputs 
%******** 
%    1  2   3   4
%   ______________
y = [C  I   R   Y];

% Imputs of State_Space Model
% ***************************
%   1      2     3                           
%  _________________
 u=[const, Gnew, dMnew ];
 
% Initial Values x0
% *****************
%    1      2      3      4      5        6        7        8       
% ___________________________________________________________________
x0=[C(1,1);I(1,1);R(1,1);Y(1,1);Y_1(1,1);R_1(1,1);R_2(1,1);R_3(1,1)];

Ns=length(t);
% Initial condition
 x_s(:,1)=x0;
% State model
 for t=2:Ns,
    x_s(:,t) = A*x_s(:,t-1) + B*u(t,:)';
  y_s(:,t-1) = CC*x_s(:,t-1);
 end
% Final values 
if t==Ns,
  y_s(:,Ns) = CC*x_s(:,Ns);
end   
t=usamacro(5:N,1);
% Defined State-Space Structure with Known Parameters from Regression
% *******************************************************************
 msr = modstruc(A,B,CC,D,K,x0);

% Packed State-Space Theta Model with Known Parameters from Regression
% ********************************************************************
 thr = ms2th(msr,'d');
 
% Residuals associeted with  state model
% ************************************** 
%... e  = resid( [ y  u ], thr);

% Simulated System (Specified in Theta Format)
% ********************************************
 yrsim_inf = predict([y u],thr);
 yrsim_1   = predict([y u],thr,1);
 yrsim_3   = predict([y u],thr,3);
 yrsim_8   = predict([y u],thr,8);
 
figure(5) 
clf
set(gcf,'Name','Compare the simulated outputs with observed ones  C, I   ( " O L S " )');
subplot(2,1,1)
plot(t,Cnew,'r.-');
hold on
plot(t,y_s(1,:)','k-')
hold on
plot(t,yrsim_inf(:,1),'g-.')
hold on
plot(t,yrsim_1(:,1),'c--')
hold on
plot(t,yrsim_3(:,1),'m-')
hold on
plot(t,yrsim_8(:,1),'y');
hold off
axis([min(t)-1,max(t)+1,min(Cnew)-90,max(Cnew)+90]);
legend('natural values','simulation y_s','prediction horizon m = inf',...
   'prediction horizon m = 1','prediction horizon m = 3',...
   'prediction horizon m = 8');
title('Real Aggrecate Persona Consumption C        ( " O L S " )');
ylabel('C_t   [ billions of 1982 dollars ]');
xlabel('Time   [ qurter ]');
grid

subplot(2,1,2)
plot(t,Inew,'r.-');
hold on
plot(t,y_s(2,:)','k-')
hold on
plot(t,yrsim_inf(:,2),'g-.')
hold on
plot(t,yrsim_1(:,2),'c--')
hold on
plot(t,yrsim_3(:,2),'m-')
hold on
plot(t,yrsim_8(:,2),'y');
hold off
axis([min(t)-1,max(t)+1,min(Inew)-40,max(Inew)+40]);
legend('natural values','simulation y_s','prediction horizon m = inf',...
   'prediction horizon m = 1','prediction horizon m = 3',...
   'prediction horizon m = 8');
title('Real Gross Domestic Investment I        ( " O L S " )');
ylabel('I_t   [ billions of 1982 dollars ]');
xlabel('Time   [ quarter ]');
grid

figure(6)
set(gcf,'Name','Compare the simulated outputs with observed ones R, Y   ( " O L S " )');
subplot(2,1,1)
plot(t,Rnew,'r.-');
hold on
plot(t,y_s(3,:)','k-')
hold on
plot(t,yrsim_inf(:,3),'g-.')
hold on
plot(t,yrsim_1(:,3),'c--')
hold on
plot(t,yrsim_3(:,3),'m-')
hold on
plot(t,yrsim_8(:,3),'y');
hold off
axis([min(t)-1,max(t)+1,min(Rnew)-1.5,max(Rnew)+1.5]);
legend('natural values','simulation y_s','prediction horizon m = inf',...
   'prediction horizon m = 1','prediction horizon m = 3',...
   'prediction horizon m = 8');
title('Interest Rate on 3-month Treasury Bills  R        ( " O L S " )');
ylabel('R_t   [ percent per year ]');
xlabel('Time  [ quarter ]');
grid

subplot(2,1,2)
plot(t,Ynew,'r.-');
hold on
plot(t,y_s(4,:)','k-')
hold on
plot(t,yrsim_inf(:,4),'g-.')
hold on
plot(t,yrsim_1(:,4),'c--')
hold on
plot(t,yrsim_3(:,4),'m-')
hold on
plot(t,yrsim_8(:,4),'y');
hold off
axis([min(t)-1,max(t)+1,min(Ynew)-90,max(Ynew)+90]);
legend('natural values','simulation y_s','prediction horizon m = inf',...
   'prediction horizon m = 1','prediction horizon m = 3',...
   'prediction horizon m = 8');
title('Real Gross Domestic Product Y        ( " O L S " )');
ylabel('Y_t   [ billions of 1982 dollars ]');
xlabel('Time   [ quarter ]');
grid

%...diary matrix1.txt % <---<<<
disp('  ');
disp('   S I M U L T A N E O U S  S T A T E - S P A C E   M O D E L')
disp('   OLS estimation on  i n d i v i d u a l   e q u a t i o n s')
disp('  ');
disp('   Matrix A')
disp(A)
disp('   Matrix B')
disp(B)
disp('   Matrix C')
disp(CC)
disp('   Matrix D')
disp(D)
disp('   Matrix K')
disp(K)
disp('  ')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%<---<<< Zde zacnete pracovat 2.10.2006 >>> 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%...diary off % <---<<<
 
% Simultaneous Estimation of Unknown Parameters of State Model
% ############################################################

% Matrix A
% 1     2     3     4     5     6     7     8     
% C_t      I_t       R_t       Y_t       Y_t-1    R_t-1     R_t-2     R_t-3
% __________________________________________________________________________
 An = [
    0         0         0         0         0         0         0         0  % 1 C_t+1
    0         0         0         0         0         0         0         0  % 2 I_t+1
    0         0         0         0         0         0         0         0  % 3 R_t+1
    0         0         0         0         0         0         0         0  % 4 Y_t+1 
    0         0         0         0         0         0         0         0  % 5 Y_t 
    0         0         0         0         0         0         0         0  % 6 R_t
    0         0         0         0         0         0         0         0  % 7 R_t-1
    0         0         0         0         0         0         0         0];% 8 R_t-2
  
% Matrix B
% 1          2        3
% const      Gnew     dMnew 
% _________________________  
 Bn = [
      0         0         0  % 1 C_t+1
      0         0         0  % 2 I_t+1
      0         0         0  % 3 R_t+1
      0         0         0  % 4 Y_t+1
      0         0         0  % 5 Y_t
      0         0         0  % 6 R_t
      0         0         0  % 7 R_t-1
      0         0         0];% 8 R_t-2

% Matrix C
%   1      2     3     4     5     6     7     8     
% C_t    I_t   R_t   Y_t  Y_t-1  R_t-1 R_t-2  R_t-3
% _________________________________________________
 Cn = [
     0     0     0     0     0     0     0     0   % 1 C_t
     0     0     0     0     0     0     0     0   % 2 I_t
     0     0     0     0     0     0     0     0   % 3 R_t
     0     0     0     0     0     0     0     0]; % 4 Y_y

% Matrix D
 Dn = [
     0     0     0
     0     0     0
     0     0     0
     0     0     0];

% Matrix K
 Kn = [
     0     0     0     0
     0     0     0     0
     0     0     0     0
     0     0     0     0
     0     0     0     0
     0     0     0     0
     0     0     0     0
     0     0     0     0];
  
% State-space structure with unknown elements
% *******************************************
 ms = [];

% Standart state-space parametrzations into the theta model format
% ****************************************************************
% Preliminary information about parameters
 thguss= [ ];  
% Preliminary information about noise    
 noisevar= [ 11.4078^2,  0         0        0;
             0           25.0965^2 0        0; 
             0           0         0.8186^2 0; 
             0           0         0         36.5044^2];   
% Theta format             
 th0 = [];
 
% Parameter Estimation of General Linear Model ( "PEM" )
% ******************************************************
 th_pem = [];;

% Transformation to state-space form
% *********************************** 
 %...[Ae,Be,Ce,De,Ke,xe0,std_Ae,std_Be,std_Ce,std_De,std_Ke,std_xe0] = []; 
 
 disp('  ');
 disp('  Simultaneous Estimated Parameters of Statate-Space Model')
 disp('  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~');
 disp('  Ae = ');
 disp(Ae)
 disp('  std_Ae =')
 disp(std_Ae)
 disp('  Be =');
 disp(Be)
 disp('  std_Be =');
 disp(std_Be);
 disp('  Ce =');
 disp(Ce);
 disp('  std_Ce =');
 disp(std_Ce)
 disp('  De =');
 disp(De)
 disp('  std_De =');
 disp(std_De);
 disp('  Ke =');
 disp(Ke);
 disp('  std_Ke =');
 disp(std_Ke)
 disp('  xe0 =')
 disp(xe0);
 disp('  std_xe0 =');
 disp(std_xe0)
 disp('  ');

figure
set(gcf,'Name','Correlations Functions of Residuals');
subplot(2,3,1)
% Residuals associeted with  state model
% ************************************** 
 e_pem  = [];

 clf
% Initial condition
 x_e(:,1)=x0;
% State model
% *********** 
 for t=2:Ns,
    x_e(:,t) = [];
  y_e(:,t-1) = [];
 end
% Final values of outputs 
if t==Ns,
  y_e(:,Ns)  = [];
end   


% Simulated System (Specified in Theta Format)
% ********************************************
 ypems_inf = predict([y u],th_pem);
 ypems_1   = predict([y u],th_pem,1);
 ypems_3   = predict([y u],th_pem,3);
 ypems_8   = predict([y u],th_pem,8);
 
t=usamacro(5:N,1);

 figure(7) 
clf
set(gcf,'Name','Compare the simulated outputs with observed ones  C, I  ( " P E M " )');
subplot(2,1,1)
plot(t,Cnew,'r.-');
hold on
plot(t,y_e(1,:)','k')
hold on
plot(t,ypems_inf(:,1),'g-.')
hold on
plot(t,ypems_1(:,1),'c--')
hold on
plot(t,ypems_3(:,1),'m-')
hold on
plot(t,ypems_8(:,1),'y');
hold off
axis([min(t)-1,max(t)+1,min(Cnew)-65,max(y_e(1,:)')+65]);
legend('natural values','simulation y_{pem}','prediction horizon m = inf',...
   'prediction horizon m = 1','prediction horizon m = 3',...
   'prediction horizon m = 8');
title('Real Aggrecate Persona Consumption C        ( " P E M " )');
ylabel('C_t   [ billions of 1982 dollars ]');
xlabel('Time   [ qurter ]');
grid

subplot(2,1,2)
plot(t,Inew,'r.-');
hold on
plot(t,y_e(2,:)','k-')
hold on
plot(t,ypems_inf(:,2),'g-.')
hold on
plot(t,ypems_1(:,2),'c--')
hold on
plot(t,ypems_3(:,2),'m-')
hold on
plot(t,ypems_8(:,2),'y');
hold off
axis([min(t)-1,max(t)+1,min(Inew)-35,max(y_e(2,:)')+35]);
legend('natural values','simulation y_{pem}','prediction horizon m = inf',...
   'prediction horizon m = 1','prediction horizon m = 3',...
   'prediction horizon m = 8');
title('Real Gross Domestic Investment I        ( " P E M " )');
ylabel('I_t   [ billions of 1982 dollars ]');
xlabel('Time   [ quarter ]');
grid

figure(8)
set(gcf,'Name','Compare the simulated outputs with observed ones R, Y ( " P E M " )');
subplot(2,1,1)
plot(t,Rnew,'r.-');
hold on
plot(t,y_e(3,:)','k-')
hold on
plot(t,ypems_inf(:,3),'g-.')
hold on
plot(t,ypems_1(:,3),'c--')
hold on
plot(t,ypems_3(:,3),'m-')
hold on
plot(t,ypems_8(:,3),'y');
hold off
axis([min(t)-1,max(t)+1,min(Rnew)-0.5,max(Rnew)+0.5]);
legend('natural values','simulation y_{pem}','prediction horizon m = inf',...
   'prediction horizon m = 1','prediction horizon m = 3',...
   'prediction horizon m = 8');
title('Interest Rate on 3-month Treasury Bills  R        ( " P E M " )');
ylabel('R_t   [ percent per year ]');
xlabel('Time  [ quarter ]');
grid

subplot(2,1,2)
plot(t,Ynew,'r.-');
hold on
plot(t,y_e(4,:)','k-')
hold on
plot(t,ypems_inf(:,4),'g-.')
hold on
plot(t,ypems_1(:,4),'c--')
hold on
plot(t,ypems_3(:,4),'m-')
hold on
plot(t,ypems_8(:,4),'y');
hold off
axis([min(t)-1,max(t)+1,min(Ynew)-90,max(y_e(4,:)')+90]);
legend('natural values','simulation y_{pem}','prediction horizon m = inf',...
   'prediction horizon m = 1','prediction horizon m = 3',...
   'prediction horizon m = 8');
title('Real Gross Domestic Product Y        ( " P E M " )');
ylabel('Y_t   [ billions of 1982 dollars ]');
xlabel('Time   [ quarter ]');
grid

ysim = y_e';

figure(9)
 set(gcf,'Name','Impulse Responses from Input  u  to Ouyput  y  ')
 impulse(th_pem,'k.')

% Computing Zeros & Poles of Theta Format Model
  zpth = th2zp(th_pem);
  
figure(10)
 set(gcf,'Name','Zeros & Poles of Model');
 pzmap(th_pem,'m')



% Estmating Frequency Response and Spectrum by Specral Analysis
  [G_s,phiv]  = spa([y u],[-1 -1]);
  [G_th,ph_t] = spa([ysim,u],[-1,-1]);

figure
 set(gcf,'Name','Frequency Response G');
 grid
 bodeplot(G_s, G_th)
 clf



