% Program ELIMOGQT
% ****************
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 REGRESNI ROVNICE
%
  y1=C;
  X1=[const C_1 Pg Pt Pop LL Y];
  [n KPRO]=size(X1);
  YPOZ=y1;
  XPOZ=X1;
  NOBD=length(y1);
  NOB1=NOBD-1; 
  NOB2=NOBD-2;
  disp( '     ')
  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 ********************
%
disp('ESTIMATOR OLS : Ordinary Least Squares Method')
% disp('ODHAD MNČ: Obyčejná Metoda Nejmenších Čtvercu')
disp('*********************************************')

% ODHAD REGRESNICH KOEFICIENTU METODOU OLS [BPAR]
  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;
   for I = 1:NOBD,
    URE2(I)=UREZ(I)*UREZ(I);
   end
  WEGHT=sqrt(URE2);

  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-KPRO);
  SIGMA=sqrt(SIGMA2);
  
% SMERODATNE ODCHYLKY PARAMETRU [SPAR] A T-STATISTIKY [TPAR]
  DIAX=diag(INVX);
  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 [DW]
   for I = 1:NOB1,
    UREZD(I)=UREZ(I+1)-UREZ(I);
   end
  UREZD2=UREZD*UREZD';
  DW=UREZD2/UREZ2;

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

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

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

% Volba určující vysvětlující proměnné dle největší t-statistiky
  for I = 1:NOBD,
  for J = 1:KPRO,
   XVAL(I,J+1)=XPOZ(I,J);
  end 
   XVAL(I,1)=YPOZ(I);
  end
% Seřazení vysvětlujících proměnných podle hodnot určující proměnné
  XYSO=sortrows(XVAL,[3]);
% [RXYS SXYS]=size(XYSO);
  
% Určení počtu vynechaných pozorování a rozdělení vzorku :
  CVYN=round(NOBD/6);
  NO=round((NOBD-CVYN)/2);
  CVYN=NOBD-2*NO;
%  
  for I = 1:NO,
  YDE1T(I)=XYSO(I,1); 
  YDE2T(I)=XYSO(NO+CVYN+I,1)';
  for J = 1:KPRO,
  XDE1(I,J)=XYSO(I,J+1);
  XDE2(I,J)=XYSO(NO+CVYN+I,J+1);
  end 
  end

  YDE1=YDE1T';
  YDE2=YDE2T';
  
% Regrese rezidua na určující proměnnou METODOU OLS v dolním vzorku
  XMAT=XDE1'*XDE1;
  YVE1=XDE1'*YDE1;
  INVX=inv(XMAT);
  BPA1=INVX*YVE1;
  YVY1=XDE1*BPA1;
  URE1=YDE1-YVY1;
  SQR1=URE1'*URE1;
  
% Regrese rezidua na určující proměnnou METODOU OLS v horním vzorku
  XMAT=XDE2'*XDE2;
  YVE2=XDE2'*YDE2;
  INVX=inv(XMAT);
  BPA2=INVX*YVE2;
  YVY2=XDE2*BPA2;
  URE2=YDE2-YVY2;  
  SQR2=URE2'*URE2;
  
% PODÍLOVÁ F-STATISTIKA A POROVNÁNÍ S KRITICKOU HODNOTOU F-ROZDĚLENÍ
  DF1=(NOBD-CVYN)/2-KPRO;
  DF2=DF1;
  RASQ=SQR2/SQR1;
  RAST=FINV(0.95,DF1,DF2);
  disp('   GOLDFELD-QUANDTova statistika')
  disp([ RASQ ])
  disp('   Tabulková hodnota F-rozdělení')
  disp([ RAST ])
  pause

  if RAST <= RASQ
  disp('   HOMOSKEDASTICITA ZAMITNUTA NA HLADINE VYZNAMNOSTI ALFA = 0,05')
  else
  disp('   HOMOSKEDASTICITA PRIJIMANA NA HLADINE VYZNAMNOSTI ALFA = 0,05')
  end 
  pause
%  return

% ******** TEST GLEJSERuv *********   
  disp(' ' )
  disp('   VYBER VYSVETLUJICI PROMENNE DOMNELE VAZANE S ABS. REZIDUY  ' )
  XX=2;
  pause
% vektor absolutních hodnot reziduí 
  AREZ=abs(UREZ);
  XTRE(:,1) = 1;
% 
% SPECIFIKACE TVARU ZAVISLOSTI
%
% LINEARNI SCHEMA
    disp(' ')
    disp('    Výsledky Glejserova testu při uvažování LINEÁRNÍ ZÁVISLOSTI :')
    for I = 1:NOBD,
    XTRE(I,2)=XPOZ(I,XX);
    end
    disp(' ')
%  outputOLST
   outputOLSX
% LOGARITMICKE SCHEMA
    disp(' ')
    disp('    Výsledky Glejserova testu při uvažování LOGARITMICKÉ ZÁVISLOSTI :')
    for I = 1:NOBD,
    XTRE(I,2)=log(XPOZ(I,XX));
    end   
    disp(' ')
%  outputOLST
   outputOLSX
% ODMOCNINNE SCHEMA
    disp(' ')
    disp('    Výsledky Glejserova testu při uvažování ODMOCNINNÉ ZÁVISLOSTI :')
    for I = 1:NOBD,
    XTRE(I,2)=sqrt(XPOZ(I,XX));
    end
    disp(' ')
%  outputOLST
   outputOLSX
  pause
  return
 
% ODHAD REGRESNICH KOEFICIENTU METODOU OLS [BPAR]

  XMAT=XTRE'*XTRE;
  YVEK=XTRE'*AREZ;
  INVX=inv(XMAT);
  BPARTR=INVX*YVEK;

% VYROVNANE HODNOTY [YVYR]
  YVYR=XTRE*BPARTR;

% ODHADNUTA REZIDUA [UREZTR] A SOUCET JEJICH CTVERCU [UREZ2TR]
  UREZTR=YPOZ-YVYR;
  UREZ2TR=UREZTR'*UREZTR;
  for I = 1:NOBD,
   URE2TR(I)=UREZTR(I)*UREZTR(I);
  end
 
% ODHAD REZIDUALNIHO ROZPTYLU [SIGMA2]
  SIGMA2=UREZ2TR/(NOBD-1);
  SIGMATR=sqrt(SIGMA2);
  
% SMERODATNE ODCHYLKY PARAMETRU [SPAR] A T-STATISTIKY [TPAR]
  DIAX=diag(INVX);
  SPARTR=SIGMATR*sqrt(DIAX);
   for I = 1:2,
    TPARTR(I)=BPARTR(I)/SPARTR(I);
   end
   
% ZOBRAZENI MEZIVYSLEDKU [parametry]
  disp('   regresni parametry')
  disp([BPARTR'])
  disp('   smerodatne odchylky parametru')
  disp([SPARTR'])
  disp('   t-statistiky parametru')
  disp([TPARTR]) 
  pause
     
 return
