/*****************************************************************************

    Function: arma.c Autoregressive function of order 1

    Description:
    Parameters: arma [-p=NN,NN,NN...] [-q=nn,nn,nn...] [-x0] [-ZSigma] [-seed] [-number]
    Return:  a list of the x1...xN values

*****************************************************************************/
#include <stdio.h>
#include <time.h>
#include <string.h>
#include <stdlib.h>
#define MAXORDER 20
void usage();

float rannorm(long *);

int main(int argc, char *argv[])
{
int p,q;
long seed;
unsigned long  i,n;
double phi[MAXORDER],xt[MAXORDER],theta[MAXORDER],zt[MAXORDER];
double x,z,m,stdev;
char *ptr;
int j;
int verbose=0;

m=0;
x=0;
p=0;
q=0;
phi[1]=0.3;
theta[1]=0;
seed=1964;
n=10;
stdev=1;

for (i = argc-1;i > 0 ;--i )
   {
/*   printf("argv[%li]:%s\n",i,argv[i]);  */
   switch (argv[i][1])
     {
      case 'm':  /* -mNNNN mean process value */
           m=atof(argv[i]+2);
           break;
      case 'x':  /* -xNNNN initial value of x */
           x=atof(argv[i]+2);
           break;
      case 'n':  /* -nNNNN number of points */
           n=atol(argv[i]+2);
           break;
      case 's':  /* -sNNNN seed value */
           seed=-atol(argv[i]+2);
           break;
      case 'v':  /* -v verbose */
           verbose=1;
           break;
      case 'p':  /* -pNNNN autoregressive coeff. list */
           p=1;
           phi[p]=strtod(argv[i]+2,&ptr);
           while(ptr[0]==',' && q < MAXORDER){
              ++ptr;
              ++p;
              phi[p]=strtod(ptr,&ptr);
           }
           break;
      case 'q':  /* -qNNNN MA coefficient list */
           q=1;
           theta[q]=strtod(argv[i]+2,&ptr);
           while(ptr[0]==',' && q < MAXORDER){
              ++ptr;
              ++q;
              theta[q]=strtod(ptr,&ptr);
           }
           break;
      case 'z':  /* -zNNNN noise stnadard deviation */
           stdev=strtod(argv[i]+2,&ptr);
           break;
      case 'h':  /* -help */
           usage();
           exit(0);
           break;
      default:
      printf("Argument %s not recognized\n",argv[i]);
      return (1);

     }
   }

 xt[0]=x;
 if (verbose){
  printf("# arma.c \n");
  printf("# dforrest@virginia.edu 1997\n");
  printf("# ARMA(%i,%i) X0=%g Zsigma=%g seed=%li n=%li\n# ",p,q,x,stdev,-seed,n);
 }
 if (p > 0)
   { for (i=1;i<=p ; ++i)
      { if(verbose)printf("phi(%li)=%g, ",i,phi[i]);
        xt[i-1]=x;
      }
   if(verbose)printf("\n# ");
   }
 if (q>0)
   { for (i=1;i<=q ;++i )
      { if(verbose)printf(" theta(%li)=%g, ",i,theta[i]);
        zt[i-1]=0;
      }
   if(verbose)printf("\n# ");
   }
 if(verbose)printf("\n");



 for (i=1;i<=n ;i++ )
    {
     z= x =  rannorm(&seed)*stdev;

     if (q>0){
       j=1;
       while (j <=q){  /* z(t-1..t-q)* theta(1..q) */
          if( i-j > 0){
               x += theta[j] * zt[(i-j)%q];
          }
          j++;
       }
       zt[i%q]= z;
     }
     if(p>0){
     j=1;
     while(j <=p ){ /* xt-(1..p)* phi(1..p) */
          if( i-j > 0){
            x += phi[j] * xt[(i-j)%p];
          }
          j++;
     }
     xt[i%p] = x;
    }
     printf("%g\n",x+m);
    }
 return 0;
}






void usage(){
     printf("arma - autoregressive moving average ARMA(p,q) process\n");
     printf("usage: arma [-x<initial value>] [-p<phi1>[,phi2]...]\n");
     printf(" -h this screen\n");
     printf(" -v verbose - print the input data\n");
     printf(" -mXXXX set the process mean\n");
     printf(" -pXXXX,XXXX,... set autoregressive phi values\n");
     printf(" -qXXXX,XXXX,... set the Moving average theta values\n");
     printf(" -sXXXXX set random number seed\n");
     printf(" -xXXXXX set the process mean\n");
     printf(" -zXXXXX set the standard deviation of the noise\n");
     printf(" -mXXXXX set the process mean\n");
     printf(" -nIIIII set the number of points\n");
     printf("    where XXXX is a floating point number\n");
     printf("    where IIII is an integer \n");
     printf("Written by David Forrest 10/18/97 - dforrest@virginia.edu\n");

     }