function fit = garch(x,order,theta0)
% fit a unit variate GARCH(p,q) model

if nargin<2 || isempty(order)
    order = [1;1];
end
if nargin<3 || isempty(theta0)
    c = var(x)*0.98;
    theta0 = [c; 0.05; zeros(order(1)-1,1);0.93;zeros(order(2)-1,1)];
end
% order = [p q]';

    Tiny = 1e-10;
    
    p = order(1);
    q = order(2);

    A = [[0 ones(1,p+q)];-eye(p+q+1)];
    b = [1-Tiny; zeros(p+q+1,1)+Tiny];
    [theta,FVAL,EXITFLAG,OUTPUT,L,G,H] =...
        fmincon(@condGARCH,theta0,A,b,[],[],[],[],[],[],x,order);
    
    fit.c1 = theta(1);
    if p > 0 
        fit.arch = theta(2:p+1);
    else
        fit.arch = [];
    end
    if q > 0
        fit.garch = theta(p+2:p+q+1);
    else
        fit.garch = [];
    end
    fit.Hessian = H;
    fit.ExitFlag = EXITFLAG;
%    fit.logl = -OUTPUT; % log likelihood

    % extract standardized residuals and volatility
    n = length(x);
    lag = max([p q]');
    sigma2 = zeros(n,1);
    epsilon = zeros(n,1);

    sigma2(1:lag) = fit.c1/(1-sum(fit.arch)-sum(fit.garch))*ones(lag,1);
    for i = lag+1:n
        tmp = fit.c1;
        if p > 0
            tmp = tmp + fit.arch'*x(i-1:-1:i-p).^2;
        end
        if q > 0
            tmp = tmp + fit.garch'*sigma2(i-1:-1:i-q);
        end
        sigma2(i) = tmp;
        epsilon(i) = x(i)/sqrt(tmp);
    end
    fit.sigma = sqrt(sigma2);
    fit.residuals = epsilon;

    function nlogl = condGARCH(theta,x,order)
    % negative conditional log likelihood of GARCH(p,q) model

    p = order(1);
    q = order(2);

    c1 = theta(1);
    if p > 0 
        arch = theta(2:p+1);
    else
        arch = [];
    end
    if q > 0
        garch = theta(p+2:p+q+1);
    else 
        garch = [];
    end
    
    n = length(x);
    lag = max([p;q]);
    xe = [zeros(lag,1);x];
    sigma2 = c1/(1-sum(arch)-sum(garch)) * ones(n+lag,1);

    for i = lag+1 : n+lag
        if p > 0
            sigma2(i) = c1 + arch'*xe(i-1:-1:i-p).^2;
        else
            sigma2(i) = c1;
        end
        if q > 0
            sigma2(i) = sigma2(i) + garch'*sigma2(i-1:-1:i-q);
        end
        if sigma2(i) <= 0
            sigma2(i) = 1e-5;
        end
    end
    sigma2 = sigma2(lag+1:n+lag);
    nlogl = - n*log(2*pi)/2 - sum(log(sigma2))/2 - sum(x.^2./(2*sigma2));
    nlogl = - nlogl;

    if isinf(nlogl) | ~isreal(nlogl) | isnan(nlogl)
        nlogl  =  1e+20;
    end

