/*****************  ARIMA (Hoeg94)  ****************************************
** ARIMA and ARIMA are procedures to simulate ARMA and ARIMA models.      **
**   Revised 12/20/94                                                     **
** The code is submitted for public non-commercial use. There are         **
** no performance guarantees. Should the code be used in a research       **
** project appropriate attribution is requested.                          **
**                                                                        **
** Written by                                                             **
**   Esben Hoeg                                                           **
**   Department of Information Science                                    **
**   The Aarhus School of Business                                        **
**   DK-8210 Aarhus V                                                     **
**   Denmark                                                              **
**                                                                        **
**   e-mail: eh@hdc.hha.dk                                                **
**                                                                        **
** Formats:                                                               **
**                                                                        **
** x = arma(T,k,variance,phi,theta,theta0,mu);                            **
**  or                                                                    **
** x = arima(T,k,variance,phi,d,theta,theta0,mu);                         **
**                                                                        **
** ar(p) is obtained from arma when theta = 0                             **
** ma(q) is obtained from arma when phi = 0                               **
**                                                                        **
** T          is number of observations                                   **
** k          is number of series (#cols in x)                            **
** variance   is the variance of the white noise                          **
** phi        is AR polynomium                                            **
** d          is the difference-parameter                                 **
**            Note that d may be fractional, so that fractionally         **
**            integrated series may be generated.                         **
** theta      is MA polynomium                                            **
** theta0     is a mean value parameter (see below)                       **
** mu         is a mean value parameter (see below)                       **
**                                                                        **
***************************************************************************/
PROC ARMA(T,k,variance,phi,theta,theta0,mu);
local help,p,q,a,x;
/***************************************************************************
**                                                                        **
** x-mu = theta0 + phi_1*(lag(x,1)-mu) + ... +phi_p*(lag(x,p)-mu)         **
**       + e + theta_1*lag(e,1) + ... + theta_q*lag(e,q);                 **
**                                                                        **
** Note that the mean of x is E(x) = theta0/(1-phi_1 - ... -phi_p) + mu   **
** (the mean may be parameterized by theta0 or mu or both)                **
**                                                                        **
**                                                                        **
** If the variance of the noise is 0, p random starting values are        **
** generated.                                                             **
**                                                                        **
** The processs becomes sinusoidal (perhaps plus noise) if                **
** the roots of phi are on the unit circle. If the roots are outside      **
** the process will die out when the variance is zero. If the roots       **
** are inside the unit circle the process will explode.                   **
***************************************************************************/
if theta==0; q=0; else; q=rows(theta); endif;
if phi==0;   p=0; else; p=rows(phi);   endif;
a     = sqrt(variance)*rndn(2*T,k);
if variance==0; a[1:p,.]=rndn(p,k); endif;
  help = lagop(a,-theta,0,1,rows(a)) + theta0;
  if p > 0;
    x = recserar(help,help[1:p,.],phi.*.ones(1,cols(help)))+mu;
  else; x = help;
  endif;
retp(trimr(x,T,0));
endp;
/**************  END of PROC ARMA *****************************************/


proc integer(d);
 /* Returns value 1 (true) if d is an integer, otherwise 0 */
 retp(d == trunc(d));
 endp;

proc delta(d);
 local help, deltalen, deltavek;
 if integer(d);
   if d<0; "integer d must be positive"; stop;
   else; deltalen = d;
   endif;
 else; deltalen = 100;  /* chosen number of weights for fractional d      */
 endif;
 help = 1-(1+d)./seqa(1,1,deltalen);
 help[1] = -help[1]; /* d */
 retp(recsercp(help,0));
endp;

PROC ARIMA(T,k,variance,phi,d,theta,theta0,mu);
local newphi, newtheta;
    if T <= rows(delta(d));
     "The number of observations is too small!"; stop;
    endif;
if d > 0;
    if phi == 0;
      newphi = delta(d);
    else;
      newphi = -rev(trimr(polymult(-rev(phi)|1,-rev(delta(d))|1),0,1));
    endif;
    retp(arma(T,k,variance,newphi,theta,theta0,mu));
elseif d < 0;
    if theta == 0;
      newtheta = -delta(-d); /* minus because of the definition of theta   */
    else;
      newtheta = rev(trimr(polymult(rev(theta)|1,-rev(delta(-d))|1),0,1));
    endif;
    retp(arma(T,k,variance,phi,newtheta,theta0,mu));
else; /* d = 0 */
    retp(arma(T,k,variance,phi,theta,theta0,mu));
endif;
endp;
/**************  END of PROC ARIMA *****************************************/


/* LAGOP is source code obtained from an earlier GAUSS shipment from Aptech*/
/* (Acronym for lag-operator)                                              */
PROC LAGOP(x,b,x0,t0,t1);
 local n, qx0, k, k0, q, s0, s1;

  q = rows(b);
 qx0 = rows(x0);

  if x0 ne 0;

    if q /= qx0;
"ERROR: If x0 is non-zero, b and x0 must have the same number
       of rows."; end;
     endif;
       k = cols(x); k0 = cols(x0);
        if k0 ne k and k0 ne 1;
"ERROR: The columns in x and x0 do not match."; end;
        endif;
       if k0 ne k; x0 = reshape(x0,k,q)'; endif;
       x = x0 | x;

     if t0 > 0;
      s0 = q + t0;
     else;
      s0 = q + 1;
     endif;

  else;

     if t0 > 0;
      s0 = t0;
     else;
      s0 = q + 1;
     endif;

 endif;

  n = rows(x);
  if t1 <= 0;
       s1 = n;
  else;
       s1 = t1;
  endif;
  if s0 > s1 or s1 > n + q - 1;
"ERROR: Illegal values specified for starting and/or ending points."; end;
  endif;

  retp( conv(1|-b, x, s0, s1) );
endp;

