From: Anthony on
Hey folks,

I'm trying to estimate an AR(4)-ARCH(3) model by minimizing the negative of the log-likelihood function.

But my "optimal" coefficients explode and are wildly huge.

I havn't included the constraint on the coefficients to ensure that the model is stationary as I am not sure how to...I think this might be the problem...

Here's what I have written:

MAIN BLOCK:
% Estimate an AR(4) model with a constant by ordinary least squares
p = 4;
[arbetas,res] = ark(Data,p);

% Engle (1982) test for autocorrelation in the residuals.
% ARCH(3) test
q = 3; % # of lags to use in ARCH test
[archstat pval archbetas] = archk(res,q);

%Estimate a AR(4)-ARCH(3) model:

theta0 = [arbetas;archbetas]; %AR and ARCH coefficients used as initial estimates in the search for optimal coefficients
Options = optimset;
thetahat = fminsearch(@likelihood,theta0,Options,Data,p,q);

**** So my thetahats explode here... ********

FUNCTIONS:
*************************************************************
function [b r] = ark(y,k)
% Purpose: estimates an AR(k) model with a constant using OLS of the form:
% y = a0 + y(t-1) a1 + y(t-2) a2 +...+ y(t-k) ak + epsilon
% Inputs: y = dependent variable vector
% k = # of lags in the model (scalar)
% Outputs: b = vector of coefficient estimates
% r = vector of residuals

% Error checking of inputs
if nargin ~= 2
error('Incorrect number of inputs');
end
if cols(y) > 1
error('y must be a vector');
end

[nobs junk] = size(y);
% Create matrix of lagged regressors of the form:
% ones + xlag(t) + xlag(t-1) + xlag(t-2) + ... xlag(t-p)
xlag = ones(nobs-k,1);
for i =1:k
xlag = [xlag y(k+1-i:end-i,1)];
end
yp = y(k+1:end); %y vector that feeds the lags
b = inv(xlag'*xlag)*xlag'*yp; % OLS coefficent estimates
r = yp - xlag*b;
end
*************************************************************
function [archstat pval b] = archk(r,k)
% Purpose: computes a test for ARCH(k) according to Engle (1982)
%
% Ho: no autoregressive conditional heterskedasticity present in the 'k'
% lagged squared series of r.

% Ha: ARCH is present in the squared series of r.
%
% Inputs: r = vector (usually regression residuals)
% k = order of ARCH to be tested (scalar)
% Outputs: archstat = ARCH(k) test statistic follows chi-squared under the
% null hypothesis
% archstat = T*R^2 ~ Chisquared with k dof
% pval = tail probability
%
% Modeled after: arch.m by Kim Baum, Dept of Economics, Boston College
% REFERENCES:
% Engle, Robert (1982), "Autoregressive Conditional Heteroskedasticity with Estimates of
% the Variance of United Kingdom Inflation", Econometrica, vol. 50, pp.
% 987-1007

% Error checking of inputs
if nargin ~= 2
error('Incorrect number of inputs');
end
if cols(r) > 1
error('r must be a vector');
end

[nobs junk] = size(r);
r2 = r.*r;
% Create matrix of lagged r^2 with a constant
r2lag = ones(nobs-k,1);
for i =1:k
r2lag = [r2lag r2(k+1-i:end-i,1)];
end
r2 = r2(k+1:end);
b = inv(r2lag'*r2lag)*r2lag'*r2; % OLS coefficent estimates
res = r2 - r2lag*b;
r2m = r2 - mean(r2);
s = res'*res; % sum of squared residuals
R2num = s; % RSS
R2denom = r2m'*r2m; % TSS
R2 = 1 - R2num/R2denom; %R-squared = 1 - RSS/TSS

archstat = length(r2)*R2;
pval = 2*(1 - chi2cdf(abs(archstat),k));
end
*************************************************************
function [LLF] = likelihood(theta,y,p,q)

%
% Purpose: Calculates the log-likelihood function of an AR(p)-ARCH(q)

% yt = rho0 + rho1*y(t-1) + rho2*y(t-2) + ... + rhop*y(t-p) + ut
% ht = alpha0 + alpha1*h(t-1) + ... + alphaq*h(t-q)

% Inputs: y = dependent variable vector
% p = # of lags of the AR portion (rho)
% q = # of lags of the ARCH portion (alpha)
% theta = vector of coefficient estimates
% theta = [rho alpha];
% Outputs: LLF = the negative of the log-likelihood function value

rho = theta(1:p+1);
alpha = theta(p+2:end);

nobs = length(y);
yp = y(p+1:end); %feeding the lags
% Matrix of lagged regressors for AR portion with constant
xlag = ones(nobs-p,1);
for i =1:p
xlag = [xlag y(p+1-i:end-i,1)];
end

res = yp - xlag*rho; %residuals
res2 = res.*res; %residuals squared, used for estimation of ARCH portion

%Build lagged squared residuals matrix for ARCH part plus a constant
lagres2 = ones(nobs-p-q+1,1); %lose p+q+1 data points because using residuals
for i=1:q
lagres2 = [lagres2, res2(q+1-i:end-i+1,1)];
end
h = lagres2*alpha;

% m = max(p,q);
% t = (p+q):nobs;
% h = h(t);
LLF = -0.5*sum(log(h(1:end-1))) - 0.5*sum((res2(q+1:end).^2)./h(1:end-1));
%-0.5ln(2pi) excluded as it includes no parameters and does not alter the
%extrema of the function.
end
*************************************************************