clc
clear all
cd 'C:\Users\Nicolás\Documents\Clases\Econometría\'


%% Heterocedasticidad

%% Especifico modelo poblacional con Heterocedasticidad
randn('seed', 123)
%Y=X*B+e

n=10000; %Tamaño muestra 1,000,000
k=5; %Número de regresores K=k+1
K=1+k; %Número de regresores incluyendo 1's

X=[ones(n,1) (10+unidrnd(100,n,k).*rand(n,k))]; %Crea vector de num aleatorios de 0 a 100 discretos
B=[10.12 12.53 5.5 2 3.2 11.23]'; %Parámetros

%Especificación de los errores:
    %Modelo multiplicativo:
    %V(ei)=ci^2 s^2 , c^2=h(z;tetha) un vector

    s=1.5; %Especifico la desviación;
    tetha=[1.2 0.02 5.02 0.42 0.12]';
    z=X(:,2:K).^2;
    c2= ones(n,1)+z*tetha; %Modelo Random Coeficient c^2= 1+z*tetha

e=randn(n,1).*(s*c2.^(0.5)); %Errores con distintas varianzas

Y=X*B+e; %Modelo especificado;

B_ols=(X'*X)^-1*X'*Y; %Parámetros estimados
e_ols=Y-X*B_ols; %Errores muestrales

e_2=e_ols.^2; %Errores al cuadrado



%% TEST GOLDFELD AND QUANT

%Idea, si los residuos son homocedásticos, la varianza de grupos debe ser
%igual.

%Requiere normalidad de los residuos.

%Paso 1: Seleccionar un regresor Xk, y ordenar {Y,X} en función de ese
%regresor Xk, el cual es fuente de heterocedasticidad (Sube Xk aumenta Volatilidad).

% Por construcción tetha=[1.2 0.02 5.02 0.42 0.12]', entonces X(:,4) es el
% escogido.

index=X(:,4);
muestra=sortrows([Y X index],(1+K+1)); %Ordeno

%Paso 2: Dividir muestra en 3 sub muestras, eliminando la central.

muestra1=muestra(1:round(n/3),:);
muestra2=muestra(round(n/3)+1:round(2*n/3),:);
muestra3=muestra(round(2*n/3)+1:n,:);

clear muestra2
%Paso 3: Estimar OLS en muestra 1 y 3.

n1=size(muestra1,1);
n3=size(muestra3,1);

y1=muestra1(:,1);
x1=muestra1(:,2:K+1);
b1=(x1'*x1)^-1*x1'*y1;

y3=muestra3(:,1);
x3=muestra3(:,2:K+1);
b3=(x3'*x3)^-1*x3'*y3;

% Paso 4: Calcular varianza de los residuos en cada submuestra.

e1=y1-x1*b1;
e3=y3-x3*b3;

v1=(e1'*e1)/(n1-K); %Varianza 1
v3=(e3'*e3)/(n3-K); %Varianza 2

% Paso 5: Calcular estadístico con v3/v1, contrastado con F(n1-K,n3-K)

%H0 v3/v1=1 (Homocedasticidad), Ha v1/v3<>1;

estadistico=v3/v1;

F=finv(0.95,n1-K,n3-K); %Punto de corte a 95%.

criterio_gq=estadistico>F; %criterio=1 -> rechazo H0, rechazo homocedasticidad.



%% TEST LM

%Paso 1: Estimar por OLS y obtener residuos e_ols
%Paso 2: Estimar la regresión: e_ols2=sigma2+d'd+u, e2=z_lm*D+u

z_lm=[ones(n,1) X(:,2:K)]; %Primera columna intercepto
e2=e_ols.^2;

%D=[sigma2 d];
D=(z_lm'*z_lm)^-1*z_lm'*e2;

%Paso 3: Calcular R2 de la regresión y Estadístico

SCR=sum((e2-z_lm*D).*(e2-z_lm*D)); % Suma cuadrática residual
SCT=sum((e2- mean(e2)).^2); %Suma cuadrados Totales
R2=1-SCR/SCT; %R cuadrado

estad=R2*n; %Estadístico

% Test de hipótesis: H0 d=0, Ha:d<>0

p=size(z_lm,2)-1; %Resto el intercepto o sigma2


xi = chi2inv(0.95,p); %Punto de corte a 95%.

criterio=estad>xi; %criterio=1 -> rechazo H0, criterio=0-> no rechazo H0




%% FGLS

%Paso 1: Estimar B_ols
%Paso 2: Estimar residuo OLS, e_ols
%Paso 3: Calcular vector de pesos w

    w=((e_ols.^2).^-1)/sum((e_ols.^2).^-1);

%Paso 4,5,6
tol=10E-10;
B_new=B_ols;
w_new=w;

for j=1:10000

B_old=B_new;
%Paso 4: Recalcular B con pesos

A1=zeros(K,K);
A2=zeros(K,1);
for i=1:n

    
    A1= w_new(i,1)*(X(i,:)'*X(i,:))+A1;
    A2= w_new(i,1)*X(i,:)'*Y(i,:)+A2;
    
end

B_new=A1^-1*A2;

%Paso 5: Recalcular w 

e_new=Y-X*B_new;
w_new=((e_new.^2).^-1)/sum((e_new.^2).^-1);

 norm(B_new-B_old,1)

%Päso 6: Seguir en el loop si no se cumple la condición 
if norm(B_new-B_old,1)<=tol 
    break
end

end

B_fgls=B_new;

%% Eicker-White

%var_gls=((X'*X)^-1*estimador*(X'*X)^-1);

%Estimador sum(e2_i*x_i'*x_i)



c=0;

for i=1:n
c=(e_ols(i,1)^2)*(X(i,:)'*X(i,:))+c;
end

estimador=c/n;


var_gls=((X'*X))^-1*estimador*((X'*X))^-1;
var_ols=(e_ols'*e_ols)/(n-K)*(X'*X)^-1;

resultado=var_ols-var_gls;

eigenvalues=eig(resultado); %Deben ser positivos ya que implica matriz semidefinida positiva.



%% 

%Y=XB+e , D^-1=P'P, D=(P'P)^-1, D=P^-1 (P')^-1, 
%PDP'=PP^-1 (P')^-1P=I
%E(ee')=sD
%E(uu')=E(Pee'P')sPDP'=sI
%PY=PXB+Pe
%B_gls=(X' D^-1 X)^-1 X' D^-1 Y

%D ?