%% Comienzo
clc, clear, % Limpia pantalla y set de variables
randn('seed', 123), %Seteo semilla para que los números aleatorios sean los mismos.
%% Temas a revisar

%1) Test de hipótesis
%1) Ley de grandes números débil
%2) Teorema central del límite univariado
%3) Teorema central del límite multivariado
%4) Test de hipótesis

%% Test de hipótesis

%Modelo y=b0+x'b1+e, quiero hacer inferencia sobre b1.

n=10000;  %Tamaño muestra
k=2; % Número de variables X, considerando los 1's
x=rand(n,1); % Genero variable x
X=[ones(n,1) x]; %Creo matriz X
e=expinv(rand(n,1)); %Creo residuos poblacionales
B=[2.1 5]'; %Parámetros
Y=X*B+e; %Modelo completo

beta=(X'*X)^-1*X'*Y; %Beta's estimados

%Test hipótesis, H0: b1=5;

%Usamos notación de restricciones lineales. Se ocupa una matriz R de
%dimensión qxk compuesta e 1 y 0, donde k es el número de variables en
%el modelo y q es el número de restricciones. r es un vector de qx1 filas.

% H0: RB=r
% Ha: RB<>r


%Ejemplo se tiene Y=b0+b1x1+b2x2+e, se quiere testear que b0=1 b1=2, por lo
%tanto se tiene, k=3, q=2. La expresión RB=r es:

%     R      B  = r
%   |1 0 0| |b0| |1|
%   |0 1 0| |b1|=|2| 
%           |b2|


R=[0 1]; %Matriz R
r=5; %Vector
q=size(r,1); %Número de restricciones
e_mco=Y-X*beta;%Errores estimados
v_mco=(e_mco'*e_mco)/(n-k);%Estimador de sigma^2


estad=(R*beta-r)'*(v_mco*R*(X'*X)^-1*R')^-1*(R*beta-r);%Estadístico
xi = chi2inv(0.95,q); %Punto de corte a 95%.
pvalue = (1-chi2cdf(estad,q)); %Pvalue

%Con 5% de significancia no rechazo H0, ya que Pvalue>0.05.


%% Ley de grandes números débil


%Secuencia iid de xi variables aleatorias con con E(xi)= u 
%y var(xi)=sigma^2, entonces media muestral converge en 
%probabilidad a la media poblacional.

n=1000000; %Tamaño
u=1.5; %Media de las variables aleatorias
s=0.9; %Desviaciones de las variables aleatorias
x=u+randn(n,1).*s;%Variables aleatorias 

%Medias muestrales para distintos tamaños de muestra

um=zeros(7,1); %vector donde guardo las medias muestrales
z=[10 100 500 1000 10000 100000 1000000]'; %vector de los tamaños de cálculo de las medias muestrales

for i=1:7
    um(i,1)=mean(x(1:z(i,1))); %Calculo medias muestrales.
end

%Calculo medias de 1 a 10, 1 a 100, 1 a 500, etc.

error=abs(um-u*ones(size(um,1),1)); % Calculo diferencias absolutas
plot(error); %Gráfico con los errores absolutos
ylabel('error');%Nombre de los ejes
xlabel('Muestra n°')

error2=(um-u*ones(size(um,1),1)).^2; % Calculo diferencias al cuadrado

plot(error2); %Gráfico con los errores cuadrados
ylabel('error');%Nombre de los ejes
xlabel('Muestra n°')



%El gráfico muestra como la diferencia entre la media muestral
%con la poblacional decrece a medida que se aumenta el tamaño 
%de la muestra.

%% Teorema central del límite univariado
clear
%Secuencia iid de xi variables aleatorias con con E(xi)= u 
%y var(xi)=sigma^2, entonces la suma de ellas
%converge en distribución a una normal.

n=1000; %Tamaño población
n2=500; %Tamaño vector z

u=1.5;
v=u^2; %2.25

z5=zeros(n2,1); %vector de variables z
z50=zeros(n2,1); %vector de variables z
z500=zeros(n2,1); %vector de variables z
z999=zeros(n2,1); %vector de variables z

for j=1:n2
    for i=1:n
x(i,j)=expinv(rand,u);
    end
end
%Cada columna es un vector de n obs.
hist(x(:,1)); 

%Calculo con distintos tamaños de muestras para ver como se acerca a una
%normal(0,1)

for j=1:n2
zz5(j,1)= 5^(0.5)*((x(1:5,j)'*ones(5,1)/5)-u)/std(x(1:5,j));
zz50(j,1)=50^(0.5)*((x(1:50,j)'*ones(50,1)/50)-u)/std(x(1:50,j));
zz500(j,1)=500^(0.5)*((x(1:500,j)'*ones(500,1)/500)-u)/std(x(1:500,j));
zz999(j,1)=999^(0.5)*((x(1:999,j)'*ones(999,1)/999)-u)/std(x(1:999,j));
end

%Cálculo de varianzas

vzz= [std(zz5)^2 std(zz50)^2 std(zz500)^2 std(zz999)^2]';

%Grafico y muestro varianzas en la parte superior.
subplot(2,2,1); hist(zz5,20);title(vzz(1))
subplot(2,2,2); hist(zz50,20);title(vzz(2))
subplot(2,2,3); hist(zz500,20);title(vzz(3))
subplot(2,2,4); hist(zz999,20);title(vzz(4))



%Si n tiende a infinito, la variable aleatoria z debería tender
%a ser una variable aleatoria distribuida normal estándar.


%% Teorema central del límite Multivariado
clear
randn('seed', 8)
%Secuencia iid de xi vectores aleatorios de tamaño k, con E(xi)= u 
%y var(xi)=S, entonces la suma de ellas
%converge en distribución a una normal.

n=1000;

%Vector de tamaño 3, defino vector de esperanzas, y
%matriz de varianzas y covarianzas.


% Se necesita generar vectores de tengan distribución normal con media 0 y
% matriz de varianzas y covarianzas sigma,

%Parte 1, crear vectores aleatorios, con cierta distribución,
%esperanza y matriz de varianza covarianza.

% x~N(0,SIGMA), ¿Cómo generar SIGMA?

%Si  x~N(0,I) -> X=w'x~N(0,w'w) donde w'w=S, S definida positiva

% Entonces tengo un vector X=~N(0,S) como dice enunciado del teorema.

%Defino primero S, y lo descompongo:

S=[4 2 1; 2 5 2; 1 2 4]; %Matriz varianza covarianza
val=eig(S); %Matriz debe ser definida positiva para poder descomponerla


w=chol(S); %Descomposición de Cholesky

w'*w % Es la matriz de varianzas y covarianzas igual a S

u=[1.5 2 6]; %Vector de esperanzas

%genero variable aleatoria con la estructura requerida,
%X=u*ones(1,n)+w'*randn(3,n); Comentado, ya que se calcula un vector en
%cada ronda.

%randn(3,n) genera un vector distribuido normal N(0,1)

%Deseo calcular z=(n^0.5)*(sum(x,2)/n-u) y ver como distribuye, calculo
%500 variables z, y estimo distintos tamaños.

z5=zeros(500,3);
z50=zeros(500,3);
z500=zeros(500,3);
zn=zeros(500,3);

%Genero vectores aleatorios, muestra de tamaño n. Necesito generar varios
%vectores ya que debo calcular la distribución de z (Medias) para distintas
%muestras.

for i=1:500
    %Nuevo Vector en cada ronda..!
 %Genero x distribuido normal e independiente.
 for j=1:n
 x=u+randn(1,3);
 X(j,:)=x*w; %3 tamaño del vector, n n° de vectores.
 end
 
  
    z5(i,:)=(5^0.5)*(sum(X(1:5,:),1)/5-u);
    z50(i,:)=(50^0.5)*(sum(X(1:50,:),1)/50-u);
    z500(i,:)=(500^0.5)*(sum(X(1:500,:),1)/500-u);
    zn(i,:)=(n^0.5)*(sum(X(1:n,:),1)/n-u);

end

%Grafico covarianzas ya que es distribución multivariada y no se puede
%graficar.

%Diferencias de las matrices de varianzas y covarianzas
%If each column is a variable, cov(X) is the covariance matrix.

dif5_n=norm(S-cov(z5),1); 
dif50_n=norm(S-cov(z50),1);
dif500_n=norm(S-cov(z500),1);
difn_n=norm(S-cov(zn),1);

%Grafico por componente

% S= Matrix 3X3 simétrica, quiero ver diferencia de diagonal y componentes
%triángulo inferior.

dif5=abs(S-cov(z5)); 
dif50=abs(S-cov(z50));
dif500=abs(S-cov(z500));
difn=abs(S-cov(zn));

%Grafico (1,1)
subplot(3,3,1);bar([dif5(1,1) dif50(1,1) dif500(1,1) dif500(1,1)]);
%Grafico (2,2)
subplot(3,3,5);bar([dif5(2,2) dif50(2,2) dif500(2,2) dif500(2,2)]);
%Grafico (3,3)
subplot(3,3,9);bar([dif5(3,3) dif50(3,3) dif500(3,3) dif500(3,3)]);

%Grafico (2,1)
subplot(3,3,4);bar([dif5(2,1) dif50(2,1) dif500(2,1) dif500(2,1)]);
%Grafico (3,1)
subplot(3,3,7);bar([dif5(3,1) dif50(3,1) dif500(3,1) dif500(3,1)]);
%Grafico (3,2)
subplot(3,3,8);bar([dif5(3,2) dif50(3,2) dif500(3,2) dif500(3,2)]);






%Calcula covarianzas

C5=cov(z5);
C50=cov(z50);
C500=cov(z500);
Cn=cov(zn);




%Calculo diferencias en varianzas columna cada vector fila componente de la
%diagonal.

pifia= zeros(3,4);
for h=1:3
        pifia(h,1) =( S(h,h)-C5(h,h))^2; 
        pifia(h,2) =( S(h,h)-C50(h,h))^2; 
        pifia(h,3) =( S(h,h)-C500(h,h))^2; 
        pifia(h,4) =( S(h,h)-Cn(h,h))^2; 
end
pifia