% E L I M O - U S A U T O
% ***********************
% Výpočet jednotlivé rovnice modelu USAuto/Maddala metodami OLS,WLS,GLS 
  close all
  clear all 
  datamodel=xlsread('dataUSAautoM');
  format 
  USAauto=datamodel;
  N=max(size(USAauto(:,:)));
%
% PROMENNE EKONOMETRICKEHO MODELU
%
%     promenne bez zpozdeni [t]
 K   =USAauto(2:N,2);  % Vzdálenost ujetá osobním autem za rok    [mile]
 C   =USAauto(2:N,3);  % Počet automobilů v USA                   [mil.ks]
 M   =USAauto(2:N,4);  % Vzdálenost ujetá autem na galon benzínu  [mile]
 Pg  =USAauto(2:N,5);  % Cena benzínu                             [index 1953]
 Pt  =USAauto(2:N,6);  % Cena veřejné dopravy                     [index 1967]
 Pop =USAauto(2:N,7);  % Obyvatelsto USA                          [mil.osob]
 L   =USAauto(2:N,8);  % Zaměstnanost USA                         [mil.osob]
 Y   =USAauto(2:N,9);  % Příjem na obyvatele USA                  [index 1958]
 t   =USAauto(2:N,1);  % Časová proměnná                          [roky]
    T=t;
%     promenna se zpozdenim 1 [t-1]
 C_1 =USAauto(1:N-1,3);  % C(t-1)
%     podílová proměnná [LL]
 LL  =L./Pop;            % podíl ekonomicky aktivního obyvatelstva
%     Transformovaná proměnná [GL]
 GPop=C.*K./(M.*Pop);    % spotřeba benzínu na obyvatele
% 
 GL  =C.*K./(M.*L);      % spotřeba benzínu na ekon.aktivního obyvatele 
%     aditivni konstanta
 const=ones(size(K)); 
%
%  SPECIFIKACE 1.REGRESNI ROVNICE
%  nobs = počet pozorování, nvar = počet vysvětlujících proměnných
%
%  y1=K;
%  X1=[const C M Pg Y];
%  [nobs nvar]=size(X1);
%  YPOZ = y1;
%  XPOZ = X1;
%  nobs=length(y1);
%  disp('  K=f(const      C        M        Pg        Y)')
%  disp('  *********************************************')
%  outputOLST
% 
% SPECIFIKACE 2.REGRESNI ROVNICE G.S.MADDALA
%
%  y1=C;
%  X1=[const C_1 Pg Pt Pop LL Y];
%  [nobs nvar]=size(X1);
%  YPOZ = y1;
%  XPOZ = X1;
%  disp(' C=f(const    C-1       Pg         Pt        Pop         LL       Y)')
%  disp(' *******************************************************************')
%  EST=1;
%  outputOLST
%  pause
% 
% SPECIFIKACE 3.REGRESNI ROVNICE G.S.MADDALA
%
%  y1=GPop;
%  X1=[const K Pg Pt Pop Y];
%  [nobs nvar]=size(X1);
%  YPOZ = y1;
%  XPOZ = X1;
%  disp(' GPop=f(const     K        Pg         Pt        Pop        Y)')
%  disp(' ************************************************************')
%  EST=1;
%  outputOLST
%  pause
%
% SPECIFIKACE 4.REGRESNI ROVNICE G.S.MADDALA
%
%  y1=GL;
%  X1=[const C_1 Pg Pt LL Y];
%  [nobs nvar]=size(X1);
%  YPOZ = y1;
%  XPOZ = X1;
%  disp(' GL=f(const     C-1       Pg         Pt        LL       Y)')
%  disp(' *********************************************************')
%  EST=1;
%  outputOLST
%  pause
% 
  disp('               SPECIFIKACE ANALYZOVANE REGRESNI ROVNICE :           ') 
  disp(' K(t)= b0 +  b1*C(t)  +  b2*M(t)  +  b3*Pg(t)  +  b4*Y(t)  +  u(t) )')
  disp(' *******************************************************************')
  y1=K;
  X1=[const C M Pg Y];
  YPOZ = y1;
  XPOZ = X1;
  [nobs nvar]=size(X1);
  
% KORELAČNÍ MATICE VYSVĚTLUJÍCÍCH PROMĚNNÝCH (BEZ VEKTORU 1)
%  
   for I = 1:nvar-1
    XCOR(:,I)=XPOZ(:,I+1);   
   end
  
   CORX=corrcoef(XCOR);  
   for I = 1:nvar-1
   for J = 1:nvar-1  
    CORR(I+1,J+1)=CORX(I,J);
   end   
   end   
   CORR(1,:)=0;
   CORR(:,1)=0;
   CORR(1,1)=1;
   disp( '     ')
   disp(' MATICE PÁROVÝCH KORELACÍ VYSVĚTLUJÍCÍCH PROMĚNNÝCH  ') 
   disp(' ***************************************************  ')
   disp([CORR])
   disp( '     ')
   pause
% 
% *********************** METODA OLS **********************
  YDEP=YPOZ; 
  XIND=XPOZ;
  disp('            ESTIMATOR OLS : Ordinary Least Squares Method')
% disp('ODHAD PMNC : Prostá metoda nejmenších čtverců')
  disp('            *********************************************')
disp('  K=f(const     C         M         Pg        Y      )')
disp('  ****************************************************')
  EST=1;
  outputWLS

  BPAROLS=BETA;
  YVYROLS=YVYR;
  UREZOLS=URTR;
  pause
  disp( '     ')
% *********************** METODA WLS **********************
  disp ('           ESTIMATOR WLS : Weighted Least Squares Method')
% disp ('ODHAD VMNC : Vážená metoda nejmenších čtverců')
  disp ('           *********************************************')
disp('  K=f(const     C         M         Pg        Y      )')
disp('  ****************************************************')
  
% STANOVENÍ RECIPR.VAH JAKO PODÍLU VYROVN.ZÁVISLE PROMĚNNÉ A ABS.REZIDUA 
   for I = 1:nobs,
%   IW(I)=abs(YVYROLS(I))./abs(UREZOLS(I));
    IW(I)=YPOZ(I);   
   end
% NORMOVÁNÍ VAH NA DIAGONÁLE MATICE V-1 NA SOUČET ROVNÝ NOBD
    CWE2=sum(IW);
    WN=IW./CWE2;
%  TRANSFORMACE PROMENNYCH
   for I = 1:nobs,
     YTR(I)=YPOZ(I)./WN(I);
   for J = 1:nvar;
     XTR(I,J)=XPOZ(I,J)./WN(I);
   end     
   end
   
   YDEP=YTR';
   XIND=XTR;
   EST=2;
   outputWLS
   
   BPARWLS=BETA;
   YVYRWLS=YVYR;
   UREZWLS=URTR;
   pause
   disp( '     ')
   
  % NORMOVÁNÍ VAH NA DIAGONÁLE MATICE V-1 NA SOUČET ROVNÝ nobs
  for I = 1:nobs,
  for J = 1:nobs,
   MAV1(I,J)=0.;  
  end   
  end
  
  for I = 1:nobs,
   IW2(I)=nobs/CWE2*IW(I)*IW(I);
   MAV1(I,I)=IW2(I);
  end     
   MAV2=MAV1;   
 
% *********************** METODA GLS se schématem AR1 ****************

  disp ('         ESTIMATOR GLS : Generalized Least Squares Method')
% disp ('ODHAD ZMNČ : Zobecněná metoda nejmenších čtverců')
  disp ('         *************************************************')
disp('  K=f(const     C         M         Pg        Y      )')
disp('  ****************************************************')
  
% VYPOCET KOEFICIENTU AUTOKORELACE 1. RADU
   UREZ=UREZOLS;
   clear UREZ1
   for I = 1:nob1,
   UREZ1(I)=UREZ(I+1);
   UREZZ(I)=UREZ(I);
   end
    
   UINV=UREZ1*UREZ1';
   UVEK=UREZ1*UREZZ';
   A1=UVEK./UINV;
%  TRANSFORMACE PROMENNYCH
    YTR(1)=YPOZ(1)*sqrt(1-A1*A1);
   for J = 1:nvar,
    XTR(1,J)=XPOZ(1,J)*sqrt(1-A1*A1);        
   end    
 %      
   for I = 2:nobs,
    YTR(I)=YPOZ(I)-A1*YPOZ(I-1);
   for J = 1:nvar;
    XTR(I,J)=XPOZ(I,J)-A1*XPOZ(I-1,J);
   end     
   end
   YDEP=YTR';
   XIND=XTR;
   EST=3;
   outputWLS
   
   BPARGLS=BETA;
   YVYRGLS=YVYR;
   UREZGLS=URTR;
   
% GRAFICKE SROVNANI REZIDUI OLS, WLS a GLS
  figure
  plot(T,UREZOLS,T,UREZWLS,T,UREZGLS)
  legend('rezidua OLS','rezidua WLS','rezidua GLS')
  title('Srovnani rezidui vypoctenych pomoci OLS, WLS a GLS')
  pause 
  return