% E L I M O - U S A U T O
% ***********************
% E L I M O - Metody WLS,GLS  USAuto
  close all
  clear all 
  datamodel=xlsread('dataUSAautoG');
  format 
  USAauto=datamodel;
  N=max(size(USAauto(:,:)));
%
% PROMENNE EKONOMETRICKEHO MODELU
%     promenne bez zpozdeni [t]
  G   =USAauto(1:N,2);  % Celková spotřeba benzínu v USA
  Pg  =USAauto(1:N,3);  % Cenový index benzínu
  Y   =USAauto(1:N,4);  % Disponibilní příjem na hlavu
  PNC =USAauto(1:N,5);  % Cenový index nových aut                
  PUC =USAauto(1:N,6);  % Cenový index ojetých aut
  PPT =USAauto(1:N,7);  % Cenový index veřejné dopravy
  PD  =USAauto(1:N,8);  % agregátní index trvanlivých komodit
  PN  =USAauto(1:N,9);  % agregátní index netrvanlivých komodit
  PS  =USAauto(1:N,10); % agregátní index služeb
  Pop =USAauto(1:N,11); % počet obyvatel v USA
  t   =USAauto(1:N,12); % pořadová časová proměnná
%
%     podílové a transformované proměnné [GP,YP,PDNS]
%
  GP  =G./Pop ;             % spotřeba benzínu na obyvatele
  YP  =GP.*Pg ;             % výdaje na benzín na obyvatele
  PDNS=(PD+PN+PS)./3 ;      % součtovyY agregátní cenový index souhrnny
% 
%    aditivni konstanta
  const=ones(size(G)); 
%
  T=t;
%
% SPECIFIKACE 1.REGRESNI ROVNICE
%
  y1=G;
  X1=[const Y  Pg PNC PUC  PPT  Pop];
  [nobs nvar]=size(X1);
  YPOZ = y1;
  XPOZ = X1;
  nobs=length(y1);
  EST=1;
% 
% SPECIFIKACE 2.REGRESNI ROVNICE
%
%  y1=GP;
%  X1=[const Y  Pg  PNC PD PN PS ];
%  [nobs nvar]=size(X1);
%  YPOZ = y1;
%  XPOZ = X1;
%  nobs=length(y1);
%  disp(' GP=f(const    Y     Pg     PNC      PD       PN       PS)')
%  disp(' *********************************************************')
%  EST=1;
% 
% SPECIFIKACE 3.REGRESNI ROVNICE
%
%  y1=GP;
%  X1=[const PUC PPT PD PN PS];
%  [nobs nvar]=size(X1);
%  YPOZ = y1;
%  XPOZ = X1;
%  nobs=length(y1);
%  disp(' YP=f(const    Y     PUC      PPT       PD       PN       PS )')
%  disp(' *********************************************************')
%  EST=1;
%
% SPECIFIKACE 4.REGRESNI ROVNICE
%
%  y1=YP;
%  X1=[const  Y  PNC  PUC  PPT  PDNS  t];
%  [nobs nvar]=size(X1);
%  YPOZ = y1;
%  XPOZ = X1;
%  nobs=length(y1);
%  disp(' YP=f(const    Y    PNC     PUC      PPT       PDNS     t)')
%  disp(' *********************************************************')%

% E L I M O - Metody WLS,GLS  USAuto
% **********************************
  disp('               SPECIFIKACE ANALYZOVANE REGRESNI ROVNICE :            ') 
  disp('  G=f(const   Y     Pg     PNC     PUC      PPT      Pop )')
  disp('  ********************************************************')
  
% 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 VMNC : Vážená metoda nejmenších čtverců')
  disp('            *********************************************')
  disp('  G=f(const   Y     Pg     PNC     PUC      PPT      Pop )')
  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('  G=f(const   Y     Pg     PNC     PUC      PPT      Pop )')
  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Ý NOBD
  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;   
%  TRA=trace(MAV1)
%  XMAT=XPOZ'*MAV1*XPOZ;
%  YVEK=XPOZ'*MAV1*YPOZ;
%  INVX=inv(XMAT);
%  BPAR=INVX*YVEK;
 
% *********************** METODA GLS **********************

  disp ('         ESTIMATOR GLS : Generalized Least Squares Method')
% disp ('ODHAD ZMNČ : Zobecněná metoda nejmenších čtverců')
  disp ('         *************************************************')
  disp('  G=f(const   Y     Pg     PNC     PUC      PPT      Pop )')
  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