1;
function [logdensity, score] = Garch(theta, data)
	mu = theta(1,:);
	omega = theta(2,:);
	alpha = abs(theta(3,:));
	beta = abs(theta(4,:));
	y = data(:,1);
	e = y - mu;
	n = rows(y);
	h = zeros(n,1);
	h(1,:) = var(y);
	for t = 2:n
	  h(t,:) = omega + alpha*e(t-1,:)^2 + beta*h(t-1,:);
	endfor
	logdensity = -log(sqrt(2*pi)) -0.5*log(h) - 0.5*(e.^2)./h;
	score = "na"; # uncomment to use numeric score vector
endfunction


load nysewk.m;
y = 100*log(close./lag(close,1));
y = y(2:rows(y),:);
plot(y); 

data = y;
names = char("constant","omega", "alpha","beta");
model = "Garch";
modelargs = "";
theta = [0.1; 0.2; 0.1; 0.9];
title = "log diff vwr, GARCH(1,1)";

ub = [0.5; 1; 1; 1];
lb = [0; 0; 0; 0];
nt = 3;
ns = 3;
rt = 0.75; # careful - this is too low for many problems
maxevals = 1e3;
neps = 5;
functol = 1e-10;
paramtol = 1e-3;
verbosity = 2; # only final results. Inc
minarg = 1;
control = { lb, ub, nt, ns, rt, maxevals, neps, functol, paramtol, verbosity,1};
theta = mle_estimate(theta, data, model, modelargs, control);
control = {-1,2};
mle_results(theta, data, model, modelargs, names, title, "");