% E L I M O - U S D A T A
% ***************
close all
load USAMACRO
format 
set(0,'DefaultLineLineWidth',2);
usamacro=USAMACRO;
N=max(size(usamacro(:,:)));
%
% PROMENNE EKONOMETRICKEHO MODELU
%
%    promenne bez zpozdeni [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)
%
%    promenne se zpozdenim 1 [t-1]
 Y_1 =usamacro(3:N-2,5);  % Y(t-1)
 R_1 =usamacro(3:N-2,4);  % R(t-1)
%
%    promenne se zpozdenim 2 [t-2]
 R_2 =usamacro(2:N-3,4);  % R(t-2)
%
%    promenne se zpozdenim 3 [t-3]
 R_3 =usamacro(1:N-4,4);  % R(t-3)
%
%    prirustkove promenne [delta t]
 dY  =Y-Y_1;              % dY=Y(t)-Y(t-1)
%
%    promenne v case 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)
%
%   diferencni promenne
 dMnew=Mnew-M;            % dM(t+1)=M(t+1)-M(t)
 dYnew=Ynew-Y;            % dY(t+1)=Y(t+1)-Y(t)
%
%    aditivni konstanta
  const=ones(size(Cnew)); % konstanta
%
  t=usamacro(5:N,1);
%
%...diary estim_us.txt %<---------<<<

% alpha is the measure of confidence (i.e. confidence  100*(1-alpha) percent) 
% alpha=0.05; 
%
% KATEGORIZACE PROMENNYCH
%
% BEZNE ENDOGENNI PROMENNE :
% C : Real Aggregate Personal CONSUMPTION
% I : Real Gross Domestic INVESTMENT
% R : Rate of Interest on 3-month Treasury bills
% Y : Real GNP
%
% PREDETERMINOVANE PROMENNE
% G : Real GOVERNMENT SPENDING
% M : Real MONEY STOCK
% Y-1 Lagged Y
% R-1 Lagged R
% R-2 Lagged R-1
% R-3 Lagged R-2
% dY  difference Y-(Y-1)
% t : time variable
%
  T=t;
%
% SPECIFIKACE 2.REGRESNI ROVNICE
%
  y1=Inew;
  X1=[const dY R-3 Ynew]; 
  [n kp]=size(X1);

% SPECIFIKACE PREDETERMINOVANYCH PROMENNYCH
  XX=[const dY R_3  dMnew (R + R_1) C];
%
  XPRD=XX;
  NOBD=length(y1);
  NOB1=NOBD-1;
  KPRO=kp;
% 
% *********************** METODA OLS ********************
%
disp('ESTIMATOR OLS : Ordinary Least Squares Method')
% disp('ODHAD MNC: obycejna Metoda Nejmensich Ctvercu')
disp('*********************************************')

% ODHAD REGRESNICH KOEFICIENTU METODOU OLS [BPAR]
  YPOZ = y1;
  XPOZ = X1;

  XMAT=XPOZ'*XPOZ;
  YVEK=XPOZ'*YPOZ;
  INVX=inv(XMAT);
  BPAR=INVX*YVEK;

% VYROVNANE HODNOTY [YVYR]
  YVYR=XPOZ*BPAR;

% ODHADNUTA REZIDUA [UREZ] A SOUCET JEJICH CTVERCU [UREZ2]
  UREZ=YPOZ-YVYR;
  UREZ2=UREZ'*UREZ;

  BPAROLS=BPAR;
  YVYROLS=YVYR;
  UREZOLS=UREZ;

% GRAF pozorovanych, vyrovnanych a rezidualnich hodnot [YPOZ,YVYR,UREZ]
  plot(T,YPOZ,T,YVYR,T,UREZ)
  legend('pozorovane','vyrovnane','rezidua')
  title('OLS : Pozorovane a vyrovnane hodnoty, rezidua rovnice')

% ODHAD REZIDUALNIHO ROZPTYLU [SIGMA2]
  SIGMA2=UREZ2/(NOBD-3);
  SIGMA=sqrt(SIGMA2);

% SMERODATNE ODCHYLKY PARAMETRU [SPAR] A T-STATISTIKY [TPAR]
  DIAX=diag(INVX);
  SPAR=SIGMA*sqrt(DIAX);
   for I = 1:3,
    TPAR(I)=BPAR(I)/SPAR(I);
   end

% ZOBRAZENI MEZIVYSLEDKU [parametry]
  disp('   regresni parametry')
  disp([BPAR'])
  disp('   smerodatne odchylky parametru')
  disp([SPAR'])
  disp('   t-statistiky parametru')
  disp([TPAR])
  pause

% KOEFICIENT DETERMINACE [R2]
  YPR=mean(YPOZ);
  YODCH=YPOZ-YPR;
  YODCH2=YODCH'*YODCH;
  R2=1.-UREZ2/YODCH2;

% DURBIN-WATSONUV KOEFICIENT
   for I = 1:NOB1,
    UREZD(I)=UREZ(I+1)-UREZ(I);
   end
  UREZD2=UREZD*UREZD';
  DW=UREZD2/UREZ2;

% von NEUMANNUV POMER
  VN=DW*NOBD/(NOBD-KPRO);

% ZOBRAZENI MEZIVYSLEDKU [vystiznost modelu]
  disp('   rezidualni rozptyl  koef.determinace')
  disp([   SIGMA           R2])

% ZOBRAZENI CHARAKTERISTIK PRIPADNE AUTOKORELACE REZIDUI
  disp('   Durbin-Watson   von Neumannuv pomer')
  disp([   DW              VN])
 pause

% *********************** METODA 2SLS ***********************

disp ('ESTIMATOR 2SLS : Two-Stage Least Squares Method')
% disp ('ODHAD 2sMNC : Dvoustupnova metoda nejmensich ctvercu')
disp ('***********************************************')

% VSTUP VSECH PREDETERMINOVANYCH PROMENNYCH (INSTRUMENTU)
  YPOZ=y1;
  XPOZ=X1;
  XPRE=XX;

% ODHAD REGRESNICH KOEFICIENTU METODOU 2SLS
 
  XMAT=XPRE'*XPRE;
  INVX=inv(XMAT);
  INMA=XPRE*INVX*XPRE';
  INMB=XPOZ'*INMA*XPOZ;
  INXX=inv(INMB);
  INXY=XPOZ'*INMA*YPOZ;
  BPAR=INXX*INXY;

% VYROVNANE HODNOTY [YVYR]
  YVYR=XPOZ*BPAR;
 
% ODHADNUTA REZIDUA A SOUCET JEJICH CTVERCU [UREZ,UREZ2]
  UREZ=YPOZ-YVYR;
  UREZ2=UREZ'*UREZ;

  BPAR2SLS=BPAR;
  YVYR2SLS=YVYR;
  UREZ2SLS=UREZ;

% GRAF pozorovanych,vyrovnanych a rezidualnich hodnot[YPOZ,YVYR,UREZ]
  plot(T,YPOZ,T,YVYR,T,UREZ)
  legend('pozorovane','vyrovnane','rezidua')
  title('2SLS : Pozorovane a vyrovnane hodnoty, rezidua rovnice')
  pause

% ODHAD REZIDUALNIHO ROZPTYLU [SIGMA2]
  SIGMA2=UREZ2/(NOBD-KPRO);
  SIGMA=sqrt(SIGMA2);

% SMERODATNE ODCHYLKY PARAMETRU [SPAR] a T-STATISTIKY [TPAR]
  DIAX=diag(INXX);
  SPAR=SIGMA*sqrt(DIAX);
   for I = 1:KPRO,
    TPAR(I)=BPAR(I)/SPAR(I);
   end

% ZOBRAZENI MEZIVYSLEDKU [parametry]
  disp('   regresni parametry')
  disp([BPAR'])
  disp('   smerodatne odchylky parametru')
  disp([SPAR'])
  disp('   t-statistiky parametru')
  disp([TPAR])
  pause

% KOEFICIENT DETERMINACE [R2]
  YPR=mean(YPOZ);
  YODCH=YPOZ-YPR;
  YODCH2=YODCH'*YODCH;
  R2=1.-UREZ2/YODCH2;

% DURBIN-WATSONUV KOEFICIENT
   for I = 1:NOB1,
    UREZD(I)=UREZ(I+1)-UREZ(I);
   end
  UREZD2=UREZD*UREZD';
  DW=UREZD2/UREZ2;

% von NEUMANNUV POMER
  VN=DW*NOBD/(NOBD-KPRO);

% ZOBRAZENI MEZIVYSLEDKU [vystiznost modelu]
  disp('   rezidualni rozptyl  koef.determinace')
  disp([   SIGMA           R2])

% ZOBRAZENI CHARAKTERISTIK PRIPADNE AUTOKORELACE REZIDUI
  disp('   Durbin-Watson   von Neumannuv pomer')
  disp([   DW              VN])
 pause
 
% GRAFICKE SROVNANI REZIDUI OLS a 2SLS
  plot(T,UREZOLS,T,UREZ2SLS)
  legend('rezidua OLS','rezidua 2SLS')
  title('Srovnani rezidui vypoctenych pomoci OLS a 2SLS')
  pause 

% *********************** METODA GIV ***********************
%
 disp ('ESTIMATOR GIV : General Instrumental Variables Method')
% disp (ODHAD OMIV : Obecna metoda instrumentalnich promennych')
 disp ('*****************************************************')

%  VSTUP VSECH INSTRUMENTALNICH PROMENNYCH
  XINS(:,1) = XX(:,1);
  XINS(:,2) = XX(:,2);
  XINS(:,3) = XX(:,4);
  XINS(:,4) = XX(:,5);
%
% XINS(:,5) = XX(:, )
% XINS(:,6) = XX(:, )

% VSTUP VSECH INSTRUMENTALNICH PROMENYCH
%
  YPOZ=y1;
  XPOZ=X1;

% ODHAD REGRESNICH KOEFICIENTU METODU IV
% 
  XMAT=XINS'*XINS;
  INVX=inv(XMAT);
  INMA=XINS*INVX*XINS';
  INMB=XPOZ'*INMA*XPOZ;
  INXX=inv(INMB);
  INXY=XPOZ'*INMA*YPOZ;
  BPAR=INXX*INXY;
%
% VYROVNANE HODNOTY [YVYR]
  YVYR=XPOZ*BPAR;
%
% ODHADNUTA REZIDUA
  UREZ=YPOZ-YVYR;
  UREZ2=UREZ'*UREZ;
  BPARIV=BPAR;
  YVYRIV=YVYR;
  UREZIV=UREZ;
%
% GRAF POZOROVANYCH, VYROVNANYCH A REZIDUALNICH HODNOT
%
  plot(T,YPOZ,T,YVYR,T,UREZ)
  legend ('pozorovane','vyrovnane','rezidua')
  title('IV : pozorovane a vyrovnane hodnoty , rezidua rovnice')
  pause

% ODHAD REZIDUALNIHO ROZPTYLU [SIGMA2]
  SIGMA2=UREZ2/(NOBD-3);
  SIGMA=sqrt(SIGMA2);

% SMERODATNE ODCHYLKY PARAMETRU [SPAR] A T-STATISTIKY [TPAR]
  DIAX=diag(INVX);
  SPAR=SIGMA*sqrt(DIAX);
   for I = 1:kp,
    TPAR(I)=BPAR(I)/SPAR(I);
   end

% ZOBRAZENI MEZIVYSLEDKU [parametry]
  disp('   regresni parametry')
  disp([BPAR'])
  disp('   smerodatne odchylky parametru')
  disp([SPAR'])
  disp('   t-statistiky parametru')
  disp([TPAR])
  pause

% KOEFICIENT DETERMINACE [R2]
  YPR=mean(YPOZ);
  YODCH=YPOZ-YPR;
  YODCH2=YODCH'*YODCH;
  R2=1.-UREZ2/YODCH2;

% DURBIN-WATSONUV KOEFICIENT
   for I = 1:NOB1,
    UREZD(I)=UREZ(I+1)-UREZ(I);
   end
  UREZD2=UREZD*UREZD';
  DW=UREZD2/UREZ2;

% von NEUMANNUV POMER
  VN=DW*NOBD/(NOBD-KPRO);

% ZOBRAZENI MEZIVYSLEDKU [vystiznost modelu]
  disp('   rezidualni rozptyl  koef.determinace')
  disp([   SIGMA           R2])

% ZOBRAZENI CHARAKTERISTIK PRIPADNE AUTOKORELACE REZIDUI
  disp('   Durbin-Watson   von Neumannuv pomer')
  disp([   DW              VN])
 pause

% GRAFICKE POROVNANI REZIDUI OLS, 2SLS a IV
 plot(T,UREZOLS,T,UREZ2SLS,T,UREZIV)
 legend('rezidua OLS', 'rezidua 2SLS','rezidua IV')
 title ('Srovnani rezidui vypoctenych pomoci OLS, 2SLS a IV')
  pause

% *********************** METODA LIML ********************************

disp('ESTIMATOR LIML : Limited Information Maximum Likelihood Method')
% disp ('ODHAD MMVOI : Metoda maximalni verohodnosti s omezenou informaci')
disp('**************************************************************')

%
% VSTUP VSECH PREDETERMINOVANYCH PROMENNYCH
%
 XPRD=XX;
 XPRM=XPRD;
 NG1=size(XPRD,2);
 XNEP(:,1)=XPRD(:,4);
 XNEP(:,2)=XPRD(:,5);
 XNEP(:,3)=XPRD(:,6);
 pause
%
% SPECIFIKACE BEZNYCH ENDOGENNICH PROMENNYCH
%
  XBEN(:,1)=y1;
  XBEN(:,2)=Ynew;
% XBEN(:,2);
% XBEN(:,3);
% XBEN(:,4);
%
 IM1=size(XBEN,2);
 IG1=3;
 IG=IG1-1;
 IM=IM1-1;
 IG2=NG1-IG1;
 IGP=IG1+1;
 KPRO=IG1;
 NOBD=length(y1);
%
 YPOZ=y1;
 XPOZ=X1;
 XPRR(:,1)=XPRD(:,1);
 XPRR(:,2)=XPRD(:,2);
 XPRR(:,3)=XPRD(:,3);
 XEND=XBEN;
 XNPR=XNEP;

% ODHAD REGRESNICH KOEFICIENTU METODOU LIML [BPAR]
 
% VYPOCET MATICE WM
  XMAT=XPRR'*XPRR;
  INVX=inv(XMAT);
  INMA=XPRR*INVX*XPRR';
  XMA1=XEND'*XEND;
  XMA2=XEND'*INMA*XEND;
  MAWM=(XMA1-XMA2)/NOBD;

% VYPOCET MATICE VM

  XMAT=XPOZ'*XPOZ;
  INVX=inv(XMAT);
  INMA=XPOZ*INVX*XPOZ';
  XMA2=XEND'*INMA*XEND;
  MAVM=(XMA1-XMA2)/NOBD;

% URCENI PRISLUSNEHO VLASTNIHO VEKTORU LVEK 
  
% SESTAVENI MATICE [(WM)^-1'*VM*(WM)^-1]
  INVX=inv(MAWM);
  XMA3=INVX'*MAVM*INVX;

% VYPOCET NEJMENSIHO VLASTNIHO CISLA LMIN MATICE [(WM)^-1'*VM*(WM)^-1]
% URCENI PRISLUSNEHO VLASTNIHO VEKTORU LVEK MATICE [VM/LMIN*WM] 
  [EIVE,EIVA]=eig(MAVM,MAWM);
%  disp('vlastni cisla a vlastni vektory') 
%  EIVE,EIVA
 
% VYPOCET SUBVEKTORU BPAR1 PARAMETRU U BEZNYCH ENDOG.PROMENNYCH
  VLVE=EIVE(:,1);
  BPA1=VLVE;
 
% NORMALIZACE SUBVEKTORU BPA1
  WNO=VLVE(1);
  BPA1(1)=-1.;
  for i = 2:IM1
   BPA1(i)=-VLVE(i)/WNO;
  end

% VYPOCET SUBVEKTORU BPA2 PARAMETRU U PREDETERMINOVANYCH PROMENNYCH
%
  XMAT=XPRR'*XPRR;
  INVX=inv(XMAT);
  INMA=INVX*XPRR';
  BPA2=-INMA*XEND*BPA1;
 
% USPORADANI DO VEKTORU BPAR
 for i = 1:IG1
  BPAR(i)=BPA2(i);
 end
 for i=1:IM1-1
  BPAR(IG1+i)=BPA1(i+1);
 end
%
% VYROVNANE HODNOTY [YVYR]
 YVYR=XPOZ*BPAR;

% ODHADNUTA REZIDUA A SOUCET JEJICH CTVRERCU [UREZ,UREZ2]
 UREZ=YPOZ-YVYR;
 UREZ2=UREZ'*UREZ;
%
 BPARLIML=BPAR;
 YVYRLIML=YVYR;
 UREZLIML=UREZ;

% GRAF pozorovanych,vyrovnanych a rezidualnich hodnot [YPOZ,YVYR,UREZ]
 plot(T,YPOZ,T,YVYR,T,UREZ)
 legend('pozorovane','vyrovnane','rezidua')
 title('Pozorovane a vyrovnane hodnoty , rezidua rovnice')
 pause

% ODHAD REZIDUALNIHO ROZPTYLU [SIGMA2]
 SIGMA2=UREZ2/(NOBD-3);
 SIGMA=sqrt(SIGMA2);

% VYPOCET ASYMPTOTICKE KOVARIANCNI MATICE ESTIMATORU LIML
%
% disp('Asymptoticka kovariancni matice estimatoru LIML')
%
% VYPOCET SUBMATICE REDUKOVANE FORMY SOUSTAVY ODPOVIDAJICI
% I-TE STRUKTURNI ROVNICI (V ENDOGENNICH PROMENNYCH)
% 
 XMAT=XPRM'*XPRM;
 INVX=inv(XMAT);
 MPRF=INVX*XPRM'*XEND;
%
% ODRIZNUTI PRVNIHO SLOUPCE MATICE PARAMETRU REDUKOVANEHO TVARU
%
 for j = 1:IM;
  M2RF(:,j)=MPRF(:,j+1);
 end
%
% ROZDELENI M2RF NA DVE MATICE [MPI1],[MPI2] PODLE PREDETERMINOVANYCH
% PROMENNYCHVE VZTAHU K I-TE ROVNICI : [PRITOMNE,NEPRITOMNE]
% 
 for i = 1:IG1;
  MPI1(i,:)=M2RF(i,:);
 end
%
 for i = 1:IG2;
  MPI2(i,:)=M2RF(IG1+i,:);
 end
%
% VYPOCET MATICE M22-M21*(M11')^-1*M21
%
 XMAT=XPRR'*XPRR;
 INVX=inv(XMAT);
 INMA=XPRR*INVX*XPRR';
%
 XMA1=XNPR'*INMA*XNPR;
 XMA2=XNPR'*XNPR;
 XMAT=XMA2-XMA1;
%
% 1. VYPOCET A DOSAZENI LEVE HORNI SUBMATICE MAAM [IMxIM]
%
 MAWX=MPI2'*XMAT*MPI2;
 DE11=inv(MAWX);
%
 for i = 1:IM;
 for j = 1:IM;
  MAAM(i,j)=DE11(i,j);
 end
 end
%
% 2. VYPOCET A DOSAZENI PRAVE HORNI SUBMATICE MAAM [IMxIG1]
%
 XMAT=XPRR'*XPRR;
 INVV=inv(XMAT);
 MAT1=INVV*XPRR'*XNPR;
%
 MADX=MPI1+MAT1*MPI2;
 DE21=MADX*DE11';
 DE12=DE21';
% 
 for i = 1:IM;
 for j = 1:IG1;
  MAAM(i,IM+j)=-DE12(i,j);
 end
 end
%
% 3.DOSAZENI LEVE DOLNI SUBMATICE MAAM [IG1xIM]
%
 for i = 1:IG1;
 for j = 1:IM;
  MAAM(IM+i,j)=-DE12(j,i);
 end
 end
%
% 4. VYPOCET A DOSAZENI PRAVE DOLNI SUBMATICE MAAM [IG1xIG1]
%
 MAWM=MADX*DE11*MADX';
 DE22=MAWM+INVX;
%
 for i = 1:IG1;
 for j = 1:IG1;
  MAAM(IM+i,IM+j)=DE22(i,j);
 end
 end
%
% VYPOCET DIAGONALNICH PRVKU MAAM = ASYMPTOTICKYCH ROZPTYLU PARAMETRU
%
 for i = 1:IM+IG1;
  SOPR(i)=SIGMA*sqrt(MAAM(i,i));
 end
%
% USPORADANI DO VEKTORU SMERODATNYCH ODCHYLEK SPAR
%
 for i = 1:IG1;
  SPAR(i)=SOPR(IM+i);
 end
 for i = 1:IM;
  SPAR(i+IG1)=SOPR(i);
 end
%
 for i = 1:IG1+IM;
  TPAR(i)=BPAR(i)/SPAR(i);
 end

% ZOBRAZENI MEZIVYSLEDKU [parametry]
  disp('   regresni parametry')
  disp([BPAR'])
  disp('   smerodatne odchylky parametru')
  disp([SPAR'])
  disp('   t-statistiky parametru')
  disp([TPAR])
  pause

% KOEFICIENT DETERMINACE [R2]
  YPR=mean(YPOZ);
  YODCH=YPOZ-YPR;
  YODCH2=YODCH'*YODCH;
  R2=1.-UREZ2/YODCH2;

% DURBIN-WATSONUV KOEFICIENT
   for I = 1:NOB1,
    UREZD(I)=UREZ(I+1)-UREZ(I);
   end
  UREZD2=UREZD*UREZD';
  DW=UREZD2/UREZ2;

% von NEUMANNUV POMER
  VN=DW*NOBD/(NOBD-KPRO);

% ZOBRAZENI MEZIVYSLEDKU [vystiznost modelu]
  disp('   rezidualni rozptyl  koef.determinace')
  disp([   SIGMA           R2])

% ZOBRAZENI CHARAKTERISTIK PRIPADNE AUTOKORELACE REZIDUI
  disp('   Durbin-Watson   von Neumannuv pomer')
  disp([   DW              VN])
 pause

% POROVNANI VYROVNANYCH HODNOT OLS,2SLS,GIV,LIML
%
disp('srovnani vyrovnanych hodnot OLS,2SLS,GIV,LIML')
disp('  OLS       2SLS      IV       LIML')
%
  YVYRV(1,:)=YVYROLS';
  YVYRV(2,:)=YVYR2SLS';
  YVYRV(3,:)=YVYRIV';
  YVYRV(4,:)=YVYRLIML';
%  YVYRV'
 pause
%
% GRAFICKE SROVNANI REZIDUI OLS, 2SLS, GIV, LIML
  plot(T,UREZOLS,T,UREZ2SLS,T,UREZIV,T,UREZLIML)
  legend('rez.OLS','rez.2SLS','rez.IV','rez.LIML')
  title('Srovnani rezidui OLS,2SLS,IV,LIML ')
%
  pause
%
