function xxx = ex01garchest
clear;
format compact;
% data load
dd = load('c:/data/ESSEX.csv');
% initial values
data = dd(:,2)/100;
irate = dd(:,3)/100/360;
data = data - irate;
mm = mean(data)
ss = std(data)
Vbar = ss^2;
b = 0.7;
c = 0.01;
par = [0. b c]; % init
garchlik(par,data,Vbar)
% constraints
parL = [-Inf 0.0001 0.0001];  % positive
parU = [ Inf 0.9999 0.9999];  % less than one
cA = [0 1 1]; % cA*par<cB
cB = [0.9999];
% 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
unitroot = par(2)+par(3)
[lik,V] = garchlik0(par, data, Vbar);
subplot(2,1,1)
plot(sqrt(V*252))
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;
if a>0
    lik = 0.; % negative likelihood
    for indx=1:N
        error = ( data(indx) - m*sqrt(VV) + 0.5*VV ) / sqrt(VV);
        lik = lik + error^2 + log(VV);
        VV = a + b*VV + c*Vbar*error^2;
        V(indx+1) = VV;
    end
else
    lik = Inf
end