%ESTIMAR MODELO GARCH
load d_ibmvwewsp6203.txt -ascii;
r = d_ibmvwewsp6203;
plot(r(:,5));
rsp = r(1:3000,5);
plot(rsp)
T = size(rsp,1);

%ESTIMACION POR REJILLA

%Estimar rsp por MCO y usar los residuos obtenidos (asegurandose de que
%no haya autocorrelacion en ellos) para estimar la varianza condicional de
%los residuos

%primero veamos la forma de la ACF y PACF de rsp
subplot(2,1,1); autocorr(rsp,30);
subplot(2,1,2); parcorr(rsp,30);
%nos sugiere un ma(1). Lo ignoramos y estimamos un AR(1)

x = [ones(2999,1) rsp(1:2999,1)];
y = rsp(2:3000,1);
beta = inv(x'*x)*(x'*y);
res = y - x*beta;
subplot(2,1,1); autocorr(res,30);
subplot(2,1,2); parcorr(res,30);

% Usar los residuos al cuadrado como proxy de ht. 

SCR = []; %poner par evitar no machacar errores
R = [];
mu = mean(rsp);
res2 = res.*res;
v = mean(res2);
LLik = [];

for i = 1:19
    for j = 1:10
        h = [];
        loglik = [];
        d1(i,1) = 0.80+i*0.010;
        d2(j,1) = 0.10+j*0.020;
        d0 = v*(1-d1(i,1)-d2(j,1));
        h(1,1) = d0 + d1(i,1)*v;
        for t = 2:T-1
          h(t,1) = d0 + d1(i,1)*h(t-1,1)+ d2(j,1)* res2(t-1,1);
          loglik = -.5*log(2*pi)-.5*log(h(t,1))-.5*res2(t,1)/h(t,1);
end
        loglik = sum(loglik);
    LLik = [LLik; loglik  d0 d1(i,1)   d2(j,1)];    
    end
end
%
a = max(LLik);
amax = max(a);
[v,p]=find(LLik==amax)
LLik(p,:)

