warning off MATLAB:divideByZero;
clear function
global T y h a ll p x;
T=length(y);
x0=[-0.2427    0.9510   -0.0307    0.4727   -5.8083    0.0311    0.0335    0.5698   -5.0257 ...
        0.0061    0.7276    5.8541];
[x,LLF,exitflag,Hess] =fminunc(@E_GARCHLLF,x0,optimset('Diagnostics','on',...
    'Display','iter','LargeScale','off','MaxFunEvals',5000,'TolFun',.0005,'HessUpdate','bfgs'),T,y,h,a);
z(1)=x(1);
z(2)=x(2);
z(3)=x(3);
z(4)=x(4);
z(5)=x(5);
z(6)=exp(x(6));
z(7)=x(7);
z(8)=x(8);
z(9)=exp(x(9));
z(10)=x(10);
z(11)=x(11);
z(12)=x(12);
p=((2^(-2/z(6)))*gamma(1/z(6))/gamma(3/z(6)))^.5;
f=p*(2^(1/z(6)))*gamma(2/z(6))/gamma(1/z(6));
h(1)=z(5);
a(1)=y(1,1)-z(7)-(z(8)*y(1,2))+(z(9)*y(1,3))+(z(10)*sin((z(11)*pi)+z(12)));
ll(1)=-(log(z(6)/p)-(.5*((abs(a(1)/(exp(h(1)/2)*p)))^z(6)))-((1+(1/z(6)))*log(2))...
    -(log(gamma(1/z(6))))-(.5*h(1)));
for j=2:T;
    h(j)=z(1)+(z(2)*h(j-1))+(z(3)*(a(j-1)/(exp(h(j-1)/2))))+(z(4)*(abs(a(j-1)/(exp(h(j-1)/2)))-f));
    a(j)=y(j,1)-z(7)-(z(8)*y(j,2))+(z(9)*y(j,3))+(z(10)*sin((z(11)*j*pi)+z(12)));
    ll(j)=-(log(z(6)/p)-(.5*((abs(a(j)/(exp(h(j)/2)*p)))^z(6)))-((1+(1/z(6)))*log(2))...
        -(log(gamma(1/z(6))))-(.5*h(j)));
end
LLF=sum(ll(:))
d=a';
c=exp(h)';
w=(a./exp(h./2));
[H,P,Qstat,CV] = lbqtest(w-mean(w), [28]', 0.05)
[H,P,ARCHstat,CV] = archtest(w-mean(w), [28]', 0.05)
for k=1:10;
    z(k)=1.00001*z(k);
p=((2^(-2/z(6)))*gamma(1/z(6))/gamma(3/z(6)))^.5;
f=p*(2^(1/z(6)))*gamma(2/z(6))/gamma(1/z(6));
h(1)=z(5);
a(1)=y(1,1)-z(7)-(z(8)*y(1,2))+(z(9)*y(1,3))+(z(10)*sin((z(11)*pi)+z(12)));
dll(1)=-(log(z(6)/p)-(.5*((abs(a(1)/(exp(h(1)/2)*p)))^z(6)))-((1+(1/z(6)))*log(2))...
    -(log(gamma(1/z(6))))-(.5*h(1)));
s(1,k)=(dll(1,k)-ll(1))/(.00001*z(k));
for j=2:T;
    h(j)=z(1)+(z(2)*h(j-1))+(z(3)*(a(j-1)/(exp(h(j-1)/2))))+(z(4)*(abs(a(j-1)/(exp(h(j-1)/2)))-f));
    a(j)=y(j,1)-z(7)-(z(8)*y(j,2))+(z(9)*y(j,3))+(z(10)*sin((z(11)*j*pi)+z(12)));
    dll(j)=-(log(z(6)/p)-(.5*((abs(a(j)/(exp(h(j)/2)*p)))^z(6)))-((1+(1/z(6)))*log(2))...
        -(log(gamma(1/z(6))))-(.5*h(j)));
    s(j,k)=(-dll(j)+ll(j))/(z(k)-(z(k)/1.00001));
end
coef(k)=z(k)/1.00001;
end
coef
se=(s'*s);
see=se\eye(10);
for j=1:10;
    stderr(j)=(see(j,j)^.5);
    tstat(j)=coef(j)/stderr(j);
    pvalue(j)=(1-tcdf(abs(tstat(j)),(1908-11)))*2;
end
pvalue
coef(11)=-.001*LLF;
coef(12)=Qstat;
coef(13)=ARCHstat;



