% test arma.m
clear all;

T = 300;
a = [0.3 0 0.2]';
b = [0 0 0 0 0 0 0 0 0 0 0 0.25]'; % MA coefficients
c = 0.01;
sigma = 0.05;

% Simulate data
n = 400;
X = zeros(n+T,1);
X(1) = 0;
p = length(a);
q = length(b);
mpq = max(p,q);
e = normrnd(0,sigma,n+T+q,1);
for i = mpq+1:T+n
    X(i) = c + e(i) + a'*X(i-1:-1:i-p) + b'*e(i-1:-1:i-q);
end

X = X(n+1:n+T);

%figure(3), plot(X)
% Estimate MA(1) Model
AR = [1 3]'; MA = [12]'; C = 1;
[theta,se,logl,res] = arma(X(1:T-10),AR,MA,C);

% Display results
theta_true = [0.7 0.2 0.25 0.01 1]';
disp('Model Estimates')
disp([theta_true theta se])

m = 10;
y = armaforecast(X(1:T-m),theta,AR,MA,C,res,m);

figure(1), plot(X,y,'.')
figure(2), plot(1:T,X,'b'), hold
plot(1:T,y,'r'), hold
