%% Generate the sample
T = 100; % Number of observations to produce
         % As you increase T, the estimated values converge to the true
         % parameter values.
omega = 0.01;
alpha = 0.07;
beta  = 0.92;
z = normrnd(0,1,[T,1]);
sigma2(1,1) = omega/(1-alpha-beta);
for i = 2:T
    sigma2(i,1) = omega + sigma2(i-1,1)*(alpha*(z(i-1,1)^2) + beta);
end
x = sqrt(sigma2).*z; % sample values
%% Initial Values for fmicnon
initialVal = [0.00001,0.03,0.96];
lb = [0;0;0];
ub = [5;1;1];
A = [0,1,1];
b = 1; % Ax <=b  [0,1,1]*[ω;α;β]<= 1
[param, fval] = fmincon(@(param) QMLE(param,x),initialVal,A,b,[],[],lb,ub,[])%,options

function LL = QMLE(param,x)
sigma2(1) = param(1)/(1 - param(2) - param(3));
LL = log(sigma2(1)) + (x(1)^2)/sigma2(1);
for i=2:length(x)
    sigma2(i) = param(1) + param(2)*x(i-1)^2 + param(3)*sigma2(i-1);
    LL = LL + log(sigma2(i)) + (x(i)^2)/sigma2(i);
end
end