function xxx = lect02garchest
clear;
format compact;
% data load
data = load('c:/dowjones.dat');
% initial values
mm = mean(data);
ss = std(data);
Vbar = ss^2;
b = 0.7;
c = 0.1;
par = [mm-.5*ss^2 b c]; % init
% constraints
parL = [-0.10 0.001 0.001];  % positive
parU = [ 0.10 0.999 0.999];  % less than one
cA = [0 1 1]; % cA*par<cB
cB = [0.999];
% optimization options
optopt = optimset('MaxIter',100,'Display','iter','LargeScale','off');
%[x, fval, exitflag, output, jacob, hess] = fminunc(@garchlik, par, optopt, data);
[x, fval, exitflag, output, lambda, jacob, hess] = fmincon(@garchlik, par, cA, cB, [], [], parL, parU, [], optopt, data, Vbar);
par = x
stderrs = diag(sqrt(inv(hess)))'
lik =-fval
[lik,V] = garchlik0(par, data, Vbar);
subplot(2,1,1)
plot(sqrt(V*12))
subplot(2,1,2)
plot(data)

function lik = garchlik(par,data,Vbar)
[lik,V] = garchlik0(par,data,Vbar);

function [lik,V] = garchlik0(par,data,Vbar)
m = par(1);
b = par(2);
c = par(3);
a = Vbar*(1-b-c);
N = size(data);
V = zeros(N+1);
VV   = Vbar;
V(1) = VV;
lik = 0.; % negative likelihood
for indx=1:N
    error = data(indx) - m + 0.5*VV;
    lik = lik + (error^2)/VV + log(VV);
    VV = a + b*VV + c*error^2;
    V(indx+1) = VV;
end