function [theta,se,logl,res] = armamle(X,theta0,p)
% estimate ARMA(p,q) model using conditional MLE
% [theta,se] = armamle(X,theta0,p)
% Input:
%   X: observations
%   theta0: initial values
%       theta0(1:p):    AR coefficients
%       theta0(p+1:p+q) MA coefficients
%       theta0(p+q+1):  standard deviation of error term
%       theta0(p+q+2):  constant term
% Output:
%   theta: estimated parameters
%   se: standard error
%   logl: log likelihood 
%   res: residuals

%[theta,logl,exitflag,output,grad,hessian] = fminunc(@lharma,theta0,[],X,p);

[theta,logl,exitflag,output,lambda,grad,hessian] = fmincon(@lharma,theta0,[],[],[],[],[],[],@stationary,[],X,p);
V = inv(hessian);
se = sqrt(diag(V));

% estimate residual
q = length(theta0)-2-p;
T = length(X);

alpha = theta(1:p);
beta = theta(p+1:p+q);
sigma = theta(1+p+q);
c = theta(2+q+p);
e = zeros(T,1);
m = max(p,q);
for ii = m+1:T
    e(ii) = X(ii) - c - alpha'*X(ii-1:-1:ii-p) ...
        - beta'*e(ii-1:-1:ii-q);
end
res = e;
logl = -logl;
function y=lharma(theta0,X,p)
%likelihood function to estimate MA(q) model

q = length(theta0)-2-p;
T = length(X);

alpha = theta0(1:p);
beta = theta0(p+1:p+q);
sigma = theta0(1+p+q);
c = theta0(2+q+p);

y = 0;
e = zeros(T,1);
m = max(p,q);
for ii = m+1:T
    e(ii) = X(ii) - c - alpha'*X(ii-1:-1:ii-p) ...
        - beta'*e(ii-1:-1:ii-q);
    y = y + e(ii)^2;
end
y = y/(2*sigma^2);
y = y + log(sigma^2)*(T-m)/2 + log(2*pi)*(T-m)/2;

if isnan(y) || isinf(y) || ~isreal(y)
    y = 1e10;
end

function [C, Ceq] = stationary(theta0,X,p)

alpha = theta0(1:p);
z = roots([-alpha(p:-1:1);1]);
C = 1-min(abs(z));
Ceq = [];



