clear
load d_ibmvwewsp6203.txt -ascii;
r = d_ibmvwewsp6203;
plot(r(:,5));
T=3000;
rsp = r(1:T,5);
plot(rsp)

%ESTIMACION POR REJILLA
SCR = [];
res = [];
mu = mean(rsp);
val = -.90:0.10:.90; % Estos son los valores alternativos de theta
% Generamos una matriz (no un vector) de residuos. 
% Cada fila son los errores de ajuste del modelo para un determinado valor de theta
n = size(val,2);
for i=1:n
    res(i,1) = rsp(1,1)-mu; 
    for j = 2:T %indice temporal
        res(i,j) = rsp(j,1)-mu +val(1,i)*res(i,j-1);
    end
end
%    Cada serie de residuos al cuadrado esta en una fila de la siguiente
%    matriz:
res_sq = res.*res;
%    Calculamos la SCR para cada  "val", sumando por filas en "res_sq" 
SCR = [SCR; val' sum(res_sq,2)];
%    Escogemos el valor de theta que minimiza la suma
[a]= find(SCR(:,2)==min(SCR(:,2)));
display('Parametros estimados mediante rejilla: mu  theta')
[mu  SCR(a,1)  sqrt(SCR(a,2)/T)]
plot(SCR(:,2));

%ESTIMACION GMM
p = sacf(rsp,12);
v = var(rsp);
mu = mean(rsp);
theta = (-1+sqrt(1-4*(p(1)^2)))/(2*p(1));
sigmae = sqrt(v/(1+theta^2));
display('Estimación:  constante  theta  DT(innovacion)')
[mu theta sigmae]

%  BOOTSTRAPPING
m = 100;
hist(rsp,m);
[freq,int] = hist(rsp,m);
fracum = [];
fracum(1,1) = freq(1,1);
for i=2:m
   fracum(i,1) = (freq(1,i)+fracum(i-1,1));
end
fracum = fracum./size(rsp,1);
plot(fracum);
%
res = [];
res(1,1) = rsp(1,1)-mu; 
for j = 2:T %indice temporal
   res(j,1) = rsp(j,1)-mu+theta*res(j-1,1);
end
sort_res=sort(res);
index = 1/T:1/T:1;
index = index';
estimaciones = [];
%  Repetir desde aqui para cambiar la estimacion
N = 20
unif = rand(T,N);
for j=1:N
draw= [];
frsp=[];
frsp(1,1) = mu;    
lagres=.0;
for i=2:T
  [v,p] = findnear(index,unif(i,j));
  %pbound = max(p,1);
  frsp(i,1)=mu+sort_res(p,1)-theta*lagres;
  lagres=sort_res(p,1);
  draw = [draw lagres];
end
%
rho = sacf(frsp,12);
vsig = var(frsp);
mu = mean(frsp);
theta = (-1+sqrt(1-4*(rho(1)^2)))/(2*rho(1));
sigmae = sqrt(vsig/(1+theta^2));
%display('Estimación:  constante  theta  DT(innovacion)')
[mu  theta  sigmae];
estimaciones = [estimaciones; mu  theta  sigmae];
end