% E L I M O - Metody WLS,GLS  USAuto
% **********************************
 close all
 clear all
 datagas=xlsread('datausagas');
 format 
 usaAuto=datagas;
 N=max(size(usaAuto(:,:)));
%
% PROMENNE EKONOMETRICKEHO MODELU
%     promenne bez zpozdeni [t]
 t   =usaAuto(2:N,1);
 K   =usaAuto(2:N,2);    % K(t)
 C   =usaAuto(2:N,3);    % C(t)
 M   =usaAuto(2:N,4);    % M(t)
 Pg  =usaAuto(2:N,5);    % Pg(t)
 Pt  =usaAuto(2:N,6);    % Pt(t)
 Pop =usaAuto(2:N,7);    % Pop(t)
 L   =usaAuto(2:N,8);    % L(t)
 Y   =usaAuto(2:N,9);    % Y(t)
 year=usaAuto(2:N,1);    % rok(t)
%     doplnené promenne
 C_1 =usaAuto(1:N-1,3);  % C(t-1)
 LL  =L./Pop;            % % zamestnanost 
%    aditivni konstanta,casova konstanta
 const=ones(size(K)); 
 T=t;
%
% SPECIFIKACE ANALYZOVANE REGRESNI ROVNICE
%
  y1=C;
  X1=[const C_1 Pg Pt Pop LL Y];
  [n kp]=size(X1);
  YPOZ = y1;
  XPOZ = X1;
  NOBD=length(y1);
  NOB1=NOBD-1;
  NOB2=NOBD-2;
  KPRO=kp;
  disp('               SPECIFIKACE ANALYZOVANE REGRESNI ROVNICE :            ') 
  disp(' C = bo*"1"  + b1*C-1  + b2*Pg   + b3*Pt  + b4*Pop   + b5*LL  +  b6*Y')
  disp(' ********************************************************************')
  
% KORELAČNÍ MATICE VYSVĚTLUJÍCÍCH PROMĚNNÝCH (BEZ VEKTORU 1)
%  
   for I = 1:KPRO-1
    XCOR(:,I)=XPOZ(:,I+1);   
   end
  
   CORX=corrcoef(XCOR);  
   for I = 1:KPRO-1
   for J = 1:KPRO-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(' C=f(const    C-1       Pg         Pt        Pop         LL       Y)')
  disp(' *******************************************************************')
  EST=1;
  outputGLS

  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(' C=f(const    C-1       Pg         Pt        Pop         LL       Y)')
  disp(' *******************************************************************')

% STANOVENÍ RECIPR.VAH JAKO PODÍLU VYROVN.ZÁVISLE PROMĚNNÉ A ABS.REZIDUA 
   for I = 1:NOBD,
     IW(I)=abs(YVYROLS(I))./abs(UREZOLS(I));
   end
   
% NORMOVÁNÍ VAH NA DIAGONÁLE MATICE V-1 NA SOUČET ROVNÝ NOBD
    CWE2=sum(IW);
    WN=IW./CWE2;
%   disp ('vahy')
%   WN
%   pause
%  TRANSFORMACE PROMENNYCH
   for I = 1:NOBD,
     YTR(I)=YPOZ(I)./WN(I);
   for J = 1:KPRO;
     XTR(I,J)=XPOZ(I,J)./WN(I);
   end     
   end
   
   YDEP=YTR';
   XIND=XTR;
   EST=2;
   outputGLS
   
   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:NOBD,
  for J = 1:NOBD,
   MAV1(I,J)=0.;  
  end   
  end
  
  for I = 1:NOBD,
   IW2(I)=NOBD/CWE2*IW(I)*IW(I);
   MAV1(I,I)=IW2(I);
  end     
   MAV2=MAV1;   

% *********************** METODA GLS **********************

  disp ('         ESTIMATOR GLS : Generalized Least Squares Method with AR(1)')
% disp ('ODHAD ZMNC : Zobecnená metoda nejmensích ctvercu' pri AR(1) )
  disp ('         ***********************************************************')
  disp(' C=f(const    C-1       Pg         Pt        Pop         LL       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:KPRO,
    XTR(1,J)=XPOZ(1,J)*sqrt(1-A1*A1);        
   end    
 %      
   for I = 2:NOBD,
    YTR(I)=YPOZ(I)-A1*YPOZ(I-1);
   for J = 1:KPRO;
    XTR(I,J)=XPOZ(I,J)-A1*XPOZ(I-1,J);
   end     
   end
   YDEP=YTR';
   XIND=XTR;
   EST=3;
   outputGLS
   
   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