display('QMLE of c+GARCH model');
T=size(y, 1);
display('insert uniform initial conditions');
coin = input('coin ');
parameters=coin*ones(4, 1);
media=parameters(1);
e=zeros(T,1);
e=y-media*ones(T,1);
parameters(2)=std(e(1:T),1)^2*(1-sum(parameters(3:4)));
iter=0;
step=0;
dLdp=zeros(4, 1);                       %--->setting fittizio
step=1;                                 %--->setting fittizio
search=norm(parameters,'fro')*ones(4,1);%--->setting fittizio



%***********************************************************
while (step*norm(search,'fro')>=10^-6*norm(parameters,'fro')),
   iter=iter+1;
   if iter>1,
      parameters=parameters+step*search;
   end,
   media=parameters(1);
   var=parameters(2:4);
   %---> reset master variables
   d=inf;
   e=zeros(T,1);
   e2=zeros(1+T,1);
   h=zeros(1+T,1);
   eps=zeros(T, 1);
   eps2=zeros(T, 1);
   hh=zeros(T, 1);
   dedm=-1;
   de2dm=zeros(1, 1+T);
   dhdm=zeros(1, 1+T);
   dhdv=zeros(3, 1+T);
   dedp=zeros(4, T);
   dhdp=zeros(4, T);
   dldp=zeros(4, T);
   %***********************************************************

   e=y-media*ones(T,1);
   h(1)=sum((e.^2))/T;
   e2(1)=sum((e.^2))/T;
   e2(2:1+T)=e.^2;

   for i=1:T,
      h(i+1)=var'*[1;e2(i);h(i)];
   end,

   likelyhood=-.5*T^-1*sum(log(h(2:1+T))+e2(2:1+T).*(h(2:1+T).^-1),1);
   L1=likelyhood;
   L2=L1;
   disp(iter),
   if iter==1, disp(likelyhood),
   else, disp([likelyhood, step, norm(dLdp, 'fro')]),
   end,

   dhdm(:, 1)=2/T*dedm*sum(e);
   de2dm(:, 1)=2/T*dedm*sum(e);
   de2dm(:, 2:1+T)=2*[dedm*e'];

   for t=1:T,
      dhdm(:, 1+t)=de2dm(:, t)*var(2)+...
         dhdm(:, t)*var(3);
      dhdv(:, 1+t)=[1;e2(t);h(t)]+...
         dhdv(:, t)*var(3);
   end;

   eps=e;
   eps2=e2(2:1+T);
   hh=h(2:1+T);
   dedp(1, :)=dedm*ones(1,T);
   dhdp=[dhdm(:, 2:1+T); dhdv(:, 2:1+T)];
   dldp=.5*repmat((eps2.*(hh.^-2)-hh.^-1)',4,1).*dhdp-...
      repmat((eps.*(hh.^-1))',4,1).*dedp;
   dLdp=T^-1*sum(dldp,2);
   %search direction****************************************
   search=inv(dldp*dldp')*sum(dldp,2);
   %********************************************************
   %in_bound conditions
   for i=2:4,
      if search(i)<0, d=min(d,-var(i-1)/search(i));
      end,
   end,

   %line search*********************************************
   lowstep=1/16;
   goldstein=0;

   if d<inf, stopstep=d;
   else, stopstep=8;
   end,


   while (lowstep<stopstep),%&(step*norm(search,'fro')>=10^-6*norm(parameters,'fro')),
      %---> reset slave variables
      par=parameters+lowstep*search;
      mediax=par(1);
      varx=par(2:4);
      ex=zeros(T,1);
      e2x=zeros(1+T,1);
      hx=zeros(1+T,1);
      epsx=zeros(T, 1);
      eps2x=zeros(T, 1);
      hhx=zeros(T, 1);
      dedmx=-1;
      de2dmx=zeros(1, T+1);
      dhdmx=zeros(1, T+1);
      dhdvx=zeros(3, T+1);
      dedpx=zeros(4, T);
      dhdpx=zeros(4, T);
      dldpx=zeros(4, T);
      dLdpx=zeros(4, 1);

      ex=y-mediax*ones(T,1);
      hx(1)=sum(ex.^2)/T;
      e2x(1)=sum(ex.^2)/T;
      e2x(2:1+T)=ex.^2;

      for i=1:T,
         hx(i+1)=varx'*[1;e2x(i);hx(i)];
      end,

      L1=-.5*T^-1*sum(log(hx(2:1+T))+e2x(2:1+T).*(hx(2:1+T).^-1),1);

      dhdmx(:, 1)=2/T*dedmx*sum(ex);
      de2dmx(:, 1)=2/T*dedmx*sum(ex);
      de2dmx(:, 2:1+T)=2*[dedmx*ex'];

      for t=1:T,
         dhdmx(:, 1+t)=de2dmx(:, t)*varx(2)+...
            dhdmx(:, t)*varx(3);
         dhdvx(:, 1+t)=[1;e2x(t);hx(t)]+...
            dhdvx(:, t)*varx(3);
      end;

      epsx=ex;
      eps2x=e2x(2:1+T);
      hhx=hx(2:1+T);
      dedpx(1, :)=dedmx*ones(1,T);
      dhdpx=[dhdmx(:, 2:1+T); dhdvx(:, 2:1+T)];
      dldpx=.5*repmat((eps2x.*(hhx.^-2)-hhx.^-1)',4,1).*dhdpx-...
         repmat((epsx.*(hhx.^-1))',4,1).*dedpx;

      dLdpx=T^-1*sum(dldpx,2);
      %******************************************************************
      if (abs(dLdpx'*search)>=.3*dLdp'*search)&(L1>=L2),
         L2=L1;
         step=lowstep;
         goldstein=1;
      end;
      lowstep=lowstep*2;

   end;

   if (goldstein==0),%&(step*norm(search,'fro')>=10^-6*norm(parameters,'fro')),
      disp('warning: step hasn''t been found correctly'),
      disp('------------------------------------------'),
      step=min(.99*d,.01);
   end;
   %************************************************************
end;

%INFORMATION MATRIX test

d2hdp=zeros(4,4,1+T);
d2hdp(1,1,1)=2;
for i=1:T,
   d2hdp(:,:,i+1)=[2*var(2),        0, de2dm(i),dhdm(i);...
                          0,        0,        0,dhdv(1,i);...
                   de2dm(i),        0,        0,dhdv(2,i);...
                    dhdm(i),dhdv(1,i),dhdv(2,i),2*dhdv(3,i)]+...
      var(3)*d2hdp(:,:,i);
end,
d2hdp=d2hdp(:,:,2:1+T);
for i=1:T,
   d2ldp(:,:,i)=1/(hh(i)^2)*(2*eps2(i)/hh(i)-1)*dhdp(:,i)*dhdp(:,i)'+...
      2/hh(i)*dedp(:,i)*dedp(:,i)'-2*eps(i)/(hh(i)^2)*(dhdp(:,i)*dedp(:,i)'+[dhdp(:,i)*dedp(:,i)']')-...
      1/hh(i)*(eps2(i)/hh(i)-1)*d2hdp(:,:,i);
end,

d2Ldp=-.5/T*sum(d2ldp,3);
for i=1:T,
   zzz(:,:,i)=dldp(:,i)*dldp(:,i)';
end,
nuova=zzz+d2ldp;

stack=zeros(10,T);
for i=1:T,
   k=1;
   for v=1:4,
      for u=1:4,
         if u<=v, stack(k,i)=nuova(u,v,i); k=k+1; end,
      end,
   end,
end,

S=(dldp*dldp')/T+d2Ldp;

k=1;
for v=1:4,
   for u=1:4,
      if u<=v, D(k,1)=S(u,v); k=k+1; end,
   end,
end,


infinitesimo=10^-6;
Dx=zeros(10,4);

for j=1:4,
   par=parameters;
   par(j)=parameters(j)+infinitesimo;
   mediax=par(1);
   varx=par(2:4);
   ex=zeros(T,1);
   e2x=zeros(1+T,1);
   hx=zeros(1+T,1);
   epsx=zeros(T, 1);
   eps2x=zeros(T, 1);
   hhx=zeros(T, 1);
   dedmx=-1;
   de2dmx=zeros(1, T+1);
   dhdmx=zeros(1, T+1);
   dhdvx=zeros(3, T+1);
   dedpx=zeros(4, T);
   dhdpx=zeros(4, T);
   dldpx=zeros(4, T);
   dLdpx=zeros(4, 1);
   d2Ldpx=zeros(4,4);
   Sx=zeros(4,4);

   ex=y-mediax*ones(T,1);
   hx(1)=sum(ex.^2)/T;
   e2x(1)=sum(ex.^2)/T;
   e2x(2:1+T)=ex.^2;

   for i=1:T,
      hx(i+1)=varx'*[1;e2x(i);hx(i)];
   end,


   dhdmx(:, 1)=2/T*dedmx*sum(ex);
   de2dmx(:, 1)=2/T*dedmx*sum(ex);
   de2dmx(:, 2:1+T)=2*[dedmx*ex'];

   for t=1:T,
      dhdmx(:, 1+t)=de2dmx(:, t)*varx(2)+...
         dhdmx(:, t)*varx(3);
      dhdvx(:, 1+t)=[1;e2x(t);hx(t)]+...
         dhdvx(:, t)*varx(3);
   end;

   epsx=ex;
   eps2x=e2x(2:1+T);
   hhx=hx(2:1+T);
   dedpx(1, :)=dedmx*ones(1,T);
   dhdpx=[dhdmx(:, 2:1+T); dhdvx(:, 2:1+T)];
   dldpx=.5*repmat((eps2x.*(hhx.^-2)-hhx.^-1)',4,1).*dhdpx-...
      repmat((epsx.*(hhx.^-1))',4,1).*dedpx;

   dLdpx=T^-1*sum(dldpx,2);

   d2hdpx=zeros(4,4,1+T);
   d2hdpx(1,1,1)=2;
   for i=1:T,
      d2hdpx(:,:,i+1)=[2*varx(2),         0, de2dmx(i),dhdmx(i);...
            0,         0,         0,dhdvx(1,i);...
            de2dmx(i),         0,         0,dhdvx(2,i);...
            dhdmx(i),dhdvx(1,i),dhdvx(2,i),2*dhdvx(3,i)]+...
         varx(3)*d2hdpx(:,:,i);
   end,
   d2hdpx=d2hdpx(:,:,2:1+T);
   for i=1:T,
      d2ldpx(:,:,i)=1/(hhx(i)^2)*(2*eps2x(i)/hhx(i)-1)*dhdpx(:,i)*dhdpx(:,i)'+...
         2/hhx(i)*dedpx(:,i)*dedpx(:,i)'-2*epsx(i)/(hhx(i)^2)*(dhdpx(:,i)*dedpx(:,i)'+[dhdpx(:,i)*dedpx(:,i)']')-...
         1/hhx(i)*(eps2x(i)/hhx(i)-1)*d2hdpx(:,:,i);
   end,
   d2Ldpx=-.5/T*sum(d2ldpx,3);
   Sx=(dldpx*dldpx')/T+d2Ldpx;
   k=1;
   for v=1:4,
      for u=1:4,
         if u<=v, Dx(k,j)=Sx(u,v); k=k+1; end,
      end,
   end,
end,
Delta=(Dx-repmat(D,1,4))./infinitesimo;
yyy=stack-((Delta*inv(d2Ldp))*dldp);
V=zeros(10,10);
for i=1:T,
   V=V+(yyy(:,i)*yyy(:,i)');
end,
V=V/T;
IM=T*D'*inv(V)*D;
IM



%clear m & n & T & coin & iter & search & step & lowstep & d & stopstep &...
%   media & var & e & e2 & h & L1 & L2 & likelyhood & goldstein & ...
%   par & mediax & varx & ex & e2x & hx & hhx & epsx & eps2x & dedmx & ...
%   de2dmx & dhdmx & dhdvx & dedpx & dhdpx & dldpx & dLdpx & i & t