function y = armaforecast(X,theta,AR,MA,C,res,n)

p = length(AR);
pmax = max(AR);
q = length(MA);
qmax = max(MA);
m = p + q + C + 1; % number of parameters.

if p>0
    alpha = theta(1:p);
    A = zeros(pmax,1);
    A(AR) = alpha;
end
if q>0
    beta = theta(p+1:p+q);
    B = zeros(qmax,1);
    B(MA) = beta;
end
if C
    c = theta(1+p+q);
else
    c = 0;
end

T = length(X);
y = [X;zeros(n,1)];
e = [res;zeros(n,1)];
%t = (1:n)';
if p>0 && q>0 
    for t = T+1:T+n
        y(t) = A'*y(t-1:-1:t-pmax)+B'*e(t-1:-1:t-qmax)+c;
    end
end
if p==0 && q>0 
    for t = T+1:T+n
        y(t) = B'*e(t-1:-1:t-qmax)+c;
    end
end
if p>0 && q==0 
    for t = T+1:T+n
        y(t) = A'*y(t-1:-1:t-pmax)+c;
    end
end


