proc (1)=newbold(ar,ma,armaser,res);
/*
newbold.src is written by Kristian Jönsson (April 10, 2003), Department of Economics, Lund University.
Contact info:  kristian.jonsson@nek.lu.se

The procedure calculates the cyclical component of a series that is to be decomposed into a permanent and a transitory part.
The calculation is performed using the methodology suggested by Newbold (1990, jME).

Inputs: 
  ar		: A px1 vector with AR coefficients.
  ma 		: A qx1 vector with MA coefficients.
  armaser	: A Tx1 vector of the stationary ARMA process.
  res		: A Tx1 residuals from the fitted ARMA series.

Output:
  cyc		: A Tx1 vector containing the cyclical component of the ARIMA(p,1,q) process.

Note:
This code can be used freely as long as proper reference is given.
No performance guarantee is made.
Bugreports are welcome.
*/
local e,A,y;
local p,q,maxpq;
local diff_pq;
local fc_first_q,rev_first_q;
local fc_to_inf;
local ct;
local i;

p=rows(ar);
q=rows(ma);
maxpq=maxc(p|q);

e=zeros(p,1);
e[1,1]=1;
A=ar'|(eye(p-1)~zeros(p-1,1));

fc_first_q=forecast(ar,ma,res,armaser,q);
fc_to_inf=zeros(rows(fc_first_q),1);

for i (maxpq,rows(armaser),1);
  rev_first_q=rev(fc_first_q[i,.]')';
  if p>q; @Extend y to include all that is needed for the forecast beyond q@
    diff_pq=p-q;
    y=(rev_first_q~armaser[i:i-diff_pq+1]')';
  elseif p<=q;
    y=fc_first_q[i,q:q-p+1]'; 
  endif;
  fc_to_inf[i,1]=e'*inv(eye(p)-A)*A*y;
endfor;

ct=sumc(fc_first_q')+fc_to_inf;

retp(ct);
endp;