program simarfit
      parameter(mt=2,ndy=16,ny=2007)
c************************************
c  Simulations of an ols fit of either an AR(p) or MA(p) model with or without
c  an additive trend or stochastic trend as described below.
c  File 3 - ARfit.out  File 2 - AR.data  File 4 - Forecast.out  File 7 - Sim.log
c========================================================
c
c  Lines for the control file SIMARFIT.CNL
c
c  The parameters & commands are found using key words whose first letter
c  may be either upper or lower case. The rest of the key word must be in
c  lower case. The symbol ! is a wildcard for the value to be set.
c  Examples: !nrun below is the wildcard for the number of runs of the program
c  using different parameters. The integer to the right of the = sign in the
c  simarfit.cnl is read by the program.
c  !sigma below is the wildcard for the standard deviation. The real number to
c  the right of = is the value of the standard deviation.
c
c  ALWAYS use 0.* rather than .* since the program searches for the number
c  to the left of the period in the parameter.
c
c  A delimter symbol * must be placed on each line.
c  Use the keywords followed by an = sign and the setting. Some settings are
c  character strings while others are real numbers.
c  IMPORTANT: Place the symbol # in a line after these parameters are set.
c  The program searches for the symbol # to delimite the parameter reads.
c  Comments on the lines of the cnl file are placed to the right of the
c  symbol *. These are not read by the program.
c
c===================================================================
c
c  Number of runs=!run (!run is an integer to be read by the program)
c
c   If run > 1 then the standard errors, their t statistics, and the
c   estimated moduli and period of the roots of the characteristic polynomial
c   of the estimated AR model are written to a file Sim.log.
c   This option is useful for production runs when run is large and will
c   only work if run > 1.
c
c============================
c
c   A character string such as 'abcd1234' are the characters enclosed by ' '.
c   These apostraphes are not typed in simarfit.cnl but are used to indicate
c   strings in the instructions. You do not have to use capital letters in
c   the control file.
c
c   Type the character string 'output data */' if you want to open a file
c   called AR.data to output the x(t) generated by the program.
c   This file will only be opened if run=1.
c
c   The data file is easily imported by Excel for making a chart (graph).
c
c   If the string 'data' is not present then this file is not opened.
c
c
c=========================
c
c  Prob threshold=!rpv (for AR t values if AR fit is used)
c   !rpv - probability level threshold for the t statistics for the AR
c    coefficients of each AR fit. If a coefficient has a t value whose
c    absolute value has a probability Pr(abs(t)) < 1 - rpv, it is deleted in
c    the next AR fit. The covariance matrix is adjusted for the subset of
c    lagged values used.
c
c   The AR model is estimated in the subroutine far which outputs the
c    po < np used as the end of at most three iterations.
c
c    If rpv=1.0 then the full AR(np) model is estimated.
c
c============================
c
c   :: Settings to select between two options for a forecasting if run=1:
c
c   Forecast='fitted model'
c   Then the fitted model is used to forecast up to 50 steps ahead.
c
c   Forecast='random walk'
c   Then the forecast is based on a random walk with a drift "trend" plus
c   an AR(p) linea dynamic model.
c
c   The random walk forecast after t = N is x(N) + c*t. It will be done only if
c   the data is detrended by first differencing.
c
c   Forecast=no
c   No forecasting is done.
c
c   If run=1 and Forecast=fitted model or Forecast=random walk then a file
c   called 'Forecast.out' is opened. This file has three columns and 50 rows.
c   The first column are the actual values of x(t) generated after the
c   sample that is generated to estimate the model parameters.
c   The second column are 50 sequential forecasts of x(t) using the true
c   AR model. The third column are 50 sequential forecasts using the
c   estimated AR parameters.
c
c============================
c
c  The loop on runs k = 1,...,nrun begins. The next set of parameter setting
c  are repeated for each run with different values.
c  Type # to delimite the end of the reads of these parameters for each run
c
c============================
c
c  Sample size=!n
c
c============================
c
c  :: Settings to control the form of the artificial data with or without an
c   additive trend.
c
c  Type the string 'ar(!p)' or 'AR(!p)' if you want to generate an AR(!p)
c  autoregressive model where !p=1 or !p is an even integer < 11

c  Type the string 'ma(!p)' or 'MA(!P)' if you want to generate an MA(!p)
c  moving average model.
c
c  IMPORTANT: The integer p must be less than or equal to 20.
c
c  The parameters of the AR(!p) process is calculated from the magnitude &
c  period of the complex roots of the model set below.
c  If an AR(1) model is selected then the lag parameter is the magnitude of
c  a real root.
c
c  The parameters for the AR model are read in the second reads loop.
c
c============================
c
c  Maximum root magnitude=!rmd
c
c   !rmd is the magnitude of the first two conjugate roots if an AR(p) model   `
c   is selected or it is the the value of the lag parameter of the AR(1).
c
c   The magnitude can GREATER than 1.0 for an MA(q) using the flipped root
c   convention where the roots have to be inside the unit circle |z|=1.
c
c   If !rmd=0 then the IMPULSE RESPONSE is calculated and is written to the
c   file AR.data.
c
c   The magnitude of the kth pair of  conjugate roots is rmd*(1-k/p)
c   for k=0,1,2,...p/2-1.
c
c   These coefficients make up the coefficients of the MA models.
c
c============================
c
c  Frequency=!fr
c
c   !fr is the frequency of the first two conjugate roots for the selected
c   AR(p) model.
c
c============================
c
c  Innovations std dev=!sig  (real)
c  !sig - standard deviation of model's stochastic input process e(t)
c
c============================
c
c  Trend type='character strings to control three options for a trend'
c
c  Trend type=none
c
c  The default is NO trend. If Trend or trend is not present in simarfit.cnl
c  then no trend is added to the AR process.
c
c  Trend type=curve
c   if you want to add a curvilinear trend to the AR process. The output is
c   x(t) = c(0) + c(1)*(tc-1/2) + c(2)*(tc**2 - tc + 1/6) + AR(1 or 2)
c   where tc = t/N.
c   The intercept c(0) and c(1) & c(2) are read from settings below.
c
c  Trend type=stochastic
c   if you want to use a stochastic trend model with p roots and a drift d.
c   Then the output is x(t) is the solution of the two equations
c
c   y(t) = d + x(t) - rho*x(t-1) & y(t) + a(1)*y(t-1) + a(p)*y(t-p) = sig*e(t)
c
c   where rho = 1 & the e(t) are unit variance innovations. This model is an
c   AR(p+1) in {x(t)} where one root is rho=1.
c   Then x(t)=d*t+y(t)+y(t-1)+...
c
c============================
c
c  If a curvilinear trend is added to the data then type
c   Trend intercept=!c
c    !c - intercept of linear trend
c
c   Linear trend slope=!s
c    !s - slope of linear part of the trend
c
c   Slope of quadratic=!q
c    !q - slope of the quadratic part of the trend
c
c  Use zeros for the intercept & slopes if you want to override the
c   switch 'curve' if it set above.
c
c  If a stochastic trend is used then type
c   Stochastic trend drift=!d
c    !d - drift constant d above
c
c  If no trend is added then leave out these command strings or put then
c  at the bottom of simarfit.cnl underneath the delimeter # or set the s
c  parameters to zero.
c
c============================
c
c  Noise model=      (for the kth run)
c   After the = sign type one of the five model types:
c   gaussian, exponential, double tailed exponential, uniform, log normal
c
c============================
c
c  :: Settings to control the method for fitting a model to the data:
c
c  Model fit='character strings to control the fitting method'
c
c  Model fit=AR(!np) or fit=ar(!np)
c   if you want to fit an AR(np) to the artificial data
c   !np - starting number of parameters of a recursive least squares AR fit
c    if np > 0. If np = 0 then this step is bypassed.
c
c  No AR fit is made if this character string is absent in the lines above #.
c
c   The residuals of the AR(np) are white if np is sufficiently large.
c   Using an AR(np) fit to generate white residuals is a linear operation
c   on the data. The AR model is estimated in the subroutine far that outputs
c   the po < np used as the end the iterations.
c
c============================
c
c  Character strings to control two options for a detrending:
c
c  Detrend='character string'
c
c  Detrend=curve
c
c  Detrend the data by fitting a linear plus quadratic curve using OLS.
c
c  Detrend='first difference'
c
c  First difference the data.
c
c  Detrend='second difference'
c
c  Second difference the data.
c
c   You can NOT use both the linear fit and differencing.
c   It is either one or the other.
c
c  Detrend='no' for no detrending or move this line to the bottom
c
c====================================
      integer run
      character*700 name
      character*80 buf
      character*50 charead
      character*20 dt,par
      character*2 delim,ot
      logical ex1
c************************************
c  Open control and summary output file
      inquire(file='simarfit.cnl',exist=ex1)
       if(.not.ex1) then
        write(6,*) 'Simarfit.cnl file does not exist in the working dire
     *ctory'
        stop
       endif
      open(1,file='simarfit.cnl',err=1,iostat=io1,status='old')
      open(3,file='ARfit.out',err=2,iostat=io2)
c  Start time
      call cpu_time(gs)
      ip=0
      do k=1,7
c  Clear buf
       do n=1,80
        buf(n:n)=' '
       enddo
c  Read line
       read(1,'(a80)') buf
       ia=index(buf,'#')
       if(ia>0) then
        exit
       endif
       ig=index(buf,'*')
       if(ig==0) then
        write(3,*) 'Place comment marker * after parameters on simarfit.
     *cnl'
        write(6,*) 'Place comment marker * after parameters on simarfit.
     *cnl'
        write(3,'(a70)')buf
        write(6,'(a70)')buf
        stop
       endif
       ig=ig-1
       name(ip+1:ip+ig)=buf(:ig)
       ip=ip+ig
      enddo
c  Number of runs
      run=1
      pv=1.
      par='runs'
      delim(1:2)='= '
      run=numread(par,name,delim,ia,ier,io)
      if(ier==1) then
       write(3,*) 'The string runs= is not in the cnl file'
       write(6,*) 'The string runs= is not in the cnl file'
       stop
      endif
c  Threshold for t value
      par='hreshold'
      pv=rnumread(par,name,delim,ia,ier,io)
      if(pv<=0. .or. pv>1.) then
       write(3,*) 't value threshold ',pv,' out of bounds'
       write(6,*) 't value threshold ',pv,' out of bounds'
       stop
      endif
      ot='nn'
      if(run==1) then
c  Check for forecast
       par='ecast'
       dt=charead(par,name,delim,ia,ier,3)
       if(ier==0) then
        i1=index(dt,'no')
        if(i1>0) then
c  No forecast
         ot(1:1)='n'
        endif
        i1=index(dt,'itt')
        if(i1>0) then
         ot(1:1)='f'
c  Forecast using the fitted model
        endif
        i1=index(dt,'walk')
        if(i1>0) then
c  Random walk forecast
         ot(1:1)='w'
        endif
       endif
       if(ot(1:1)/='n') then
c  Open file Forecast.out file
        open(4,file='Forecast.out',err=3,iostat=io3)
       endif
        elseif(run>1) then
         ot(2:2)='r'
         open(7,file='Sim.log',err=4,iostat=io4)
      endif
      is=index(name,'data')
      if(is>0 .and. run==1) then
c  Data file will be opened in loop
       open(2,file='AR.data',err=5,iostat=io5)
       ot(2:2)='d'
      endif
c-------
c  Start loop on runs
      do nr=1,run
       if(run>1) then
        write(3,'(15x,''|'',14(''=''),''|'')')
        write(3,'(15x,''| Run No. '',i3,''  |'')') nr
        write(3,'(15x,''|'',14(''=''),''|''/)')
       endif
       call loop(nr,run,pv,ot,3)
       call flush(3)
       if(ot(2:2)=='r') then
        call flush(7)
       endif
      enddo
c+++++++++++++++++++++++++++++++++++++++++++
c  Time and date of run & run time
      call head(mt,ndy,ny,gs,rt,3)
c+++++++++++++++++++++++++++++++++++++++++++
      stop
  1   write(6,*)'**** Error on opening control file 1 ',io1
      stop
  2   write(6,*)'**** Error on opening output file ',io2
      stop
  3   write(6,*)'**** Error on opening Forecast.out file ',io3
      stop
  4   write(6,*)'**** Error on opening Sim.log file ',io4
      stop
  5   write(6,*)'**** Error on opening file ar.data ',io5
      stop
      end
c
c===========================================
c
      subroutine loop(nr,run,pv,ot,io)
c********************************************************
c  Open data file and read data
c========================================================
      allocatable::a,x,s,e,r,c,v,in,mag,freq,root
      real x(:),s(:),e(:),r(:),a(:),ts(0:2),d(0:2)
      real*8 c(:),v(:),mag(:),freq(:),pi2,s1
      complex*16 root(:)
      integer in(:),p,q,po,run,seed,t
      character*700 name
      character*80 buf
      character*50 charead
      character*25 model
      character*20 dt,par
      character*2 delim,ot
      character ch,dtrend,init,mod
      intent(in) nr,run,pv,ot,io
c********************************************
      ip=0
      do k=1,17
c  Clear buf
       do n=1,80
        buf(n:n)=' '
       enddo
c  Read line
       read(1,'(a80)') buf
       ia=index(buf,'#')
       if(ia>0) then
        exit
       endif
       ig=index(buf,'*')
       if(ig==0) then
        write(io,*) 'Place comment marker * after parameters on simarfit
     *.cnl'
        write(io,'(a70)') buf
        stop
       endif
       ig=ig-1
       name(ip+1:ip+ig)=buf(:ig)
       ip=ip+ig
      enddo
c  Sample size
      par='size'
      delim(1:2)='= '
      n=numread(par,name,delim,ia,ier,io)
      if(ier>0) then
       write(io,*)'Error in reading sample size'
       write(6,*)'Error in reading sample size'
       stop
      endif
c  Model type
      par='model'
      dt=charead(par,name,delim,ia,ier,3)
      if(ier>0) then
       write(io,*)'Error in reading model type'
       write(6,*)'Error in reading model type'
       stop
      endif
      if(ia>0) then
        ig=index(name,'gauss')
        id=index(name,'aile')
        ie=index(name,'exp')
        iu=index(name,'unif')
        iln=index(name,'log')
        if(ig>0) then
         ch='g'
         model='gaussian'
          elseif(id>0) then
           ch='d'
            model='double tailed exponential'
             elseif(ie>0) then
              ch='e'
              model='exponential'
               elseif(iu>0) then
              ch='u'
              model='uniform'
               elseif(iln>0) then
                ch='l'
                model='log normal'
                 else
                  write(io,*) 'Model control letter is incorrect ',dt
                  stop
        endif
         elseif(ia==0) then
          write(io,'(/7x,''Error in read of innovations model type'')')
          write(6,'(/7x,''Error in read of innovations model type'')')
         stop
      endif
c  p - AR lag, q - MA lag
      p=0
      q=0
      par='lag'
      ar=numread(par,name,delim,ia,ier,io)
      if(ier>0) then
       write(io,'(/''Error in reading lag value for AR Or MA model''/)')
       write(6,'(/''Error in reading lag value for AR Or MA model''/)')
       stop
      endif
      i1=index(name(ia-11:ia-5),'AR')
      i2=index(name(ia-11:ia-5),'ar')
      i3=index(name(ia-11:ia-5),'MA')
      i4=index(name(ia-11:ia-5),'ma')
      if(i1>0.or.i2>0) then
       p=ar
       elseif(i3>0.or.i4>0) then
        q=ar
         else
          write(io,*)'The characters AR or ar or MA or ma are missing in
     * the Simarfit.cnl'
          write(6,*)'The characters AR or ar or MA or ma are missing in
     *the Simarfit.cnl'
          stop
      endif
c  Trap p or q
      if(p>20) then
       write(io,'(/2x,''Lag p for AR(p) > 20 '',i3)')p
       write(6,'(/2x,''Lag p for AR(p) > 20 '',i3)')p
       stop
      endif
      if(q>20) then
       write(io,'(/2x,''Lag q for MA(q) > 20 '',i3)')q
       write(6,'(/2x,''Lag q for MA(q) > 20 '',i3)')q
       stop
      endif
      delim(1:2)='= '
c  Trend model type
      par='type'
      dt=charead(par,name,delim,ia,ier,3)
      if(ier==0) then
       i1=index(dt,'no')
       if(i1>0) then
        mod='n'
       endif
       i1=index(dt,'urve')
       if(i1>0) then
c  Curvilinear trend
        mod='t'
       endif
       i1=index(dt,'stoch')
       if(i1>0) then
c  Stochastic trend
        mod='s'
       endif
        else
         write(io,'(/7x,''Error in reading the trend model in simarfit.c
     *nl'')')
         write(6,'(/7x,''Error in reading the trend model in simarfit.cn
     *l'')')
         stop
       if(mod/='n'.and.mod/='t'.and.mod/='s') then
        write(io,'(/7x,''Wrong input for the type of trend model''/)')
        write(io,'(7x,''Use none or curve or stochastic in simarfit.cnl'
     *')')
        write(6,'(/7x,''Wrong input for the type of trend model''/)')
        write(6,'(7x,''Use none or curve or stochastic in simarfit.cnl''
     *)')
        stop
       endif
      endif
c  Trap p=0 for stochastic trend
      if(mod=='s'.and.p==0) then
       write(io,'(/2x,''Set a positive lag for the AR(p) in the stochast
     *tic trend p = '',i1)')p
       write(6,'(/2x,''Set a positive lag for the AR(p) in the stochasti
     *ic trend p = '',i1)')p
       stop
      endif
c  AR fit max lag
      np=0
      par='fit'
      dt=charead(par,name,delim,ia,ier,3)
      if(ier==0) then
c  AR fit of generated data
       par='R'
       delim(1:2)='()'
       np=numread(par,dt,delim,ia,ier,3)
       if(ier>0) then
        par='r'
        np=numread(par,dt,delim,ia,ier,3)
       endif
        else
         np=0
      endif
      if(np>100) then
       write(io,'(/i3,'' > 100. Use a smaller lag value'')')np
       write(6,'(/i3,'' > 100. Use a smaller lag value'')')np
       stop
      endif
c  Root maximum magnitude
      delim(1:2)='= '
      par='nitude'
      rmd=rnumread(par,name,delim,ia,ier,io)
      if(ier>0) then
       write(io,*)'Error in reading root magnitude'
       write(6,*)'Error in reading root magnitude'
       stop
      endif
      if(rmd<=0.0) then
       write(6,'(/2x,''Max magnitude of the root '',f7.3,'' <= 0'')')rmd
       write(3,'(/2x,''Max magnitude of the root '',f7.3,'' <= 0'')')rmd
       stop
      endif
      if(rmd>1.0.and.p>0) then
       write(6,'(/2x,''Magnitude of 1st AR root '',f7.3,'' > 1.0'')')rmd
       write(3,'(/2x,''Magnitude of 1st AR root '',f7.3,'' > 1.0'')')rmd
       stop
      endif
c  Frequency of the first root
      par='quency'
      fr=rnumread(par,name,delim,ia,ier,io)
      if(ier>0) then
       write(io,*)'Error in reading root frequency'
       write(6,*)'Error in reading root frequency'
       stop
      endif
      if(fr<=0.0) then
       write(6,'(/2x,''Frequency of 1st root '',f7.3,'' <= 0'')')fr
       write(io,'(/2x,''Frequency of 1st root '',f7.3,'' <= 0'')')fr
       stop
      endif
      if(fr>0.5) then
       write(6,'(/2x,''Frequency of 1st root '',f7.3,'' > 0.5'')')fr
       write(io,'(/2x,''Frequency of 1st root '',f7.3,'' > 0.5'')')fr
       stop
      endif
c  Standard deviation for e(t)
      par='dev'
      sig=rnumread(par,name,delim,ia,ier,io)
      if(sig<0.) then
       write(io,*) 'sig for e(t)',sig,' < 0'
       write(6,*) 'sig for e(t)',sig,' < 0'
       stop
      endif
c  Detrending method
      dtrend='n'
      ifd=0
      par='etrend'
      dt=charead(par,name,delim,ia,ier,3)
c  Switches for least squares fit of trend, first or second difference of data
      if(ier==0) then
       i1=index(dt,'urve')
       if(i1>0) then
c  No differencing
        dtrend='n'
       endif
       i1=index(dt,'urve')
       if(i1>0) then
        dtrend='d'
       endif
       i1=index(dt,'rst')
       if(i1>0) then
c  1st difference
        ifd=1
       endif
       i1=index(dt,'ond')
       if(i1>0) then
c  2nd difference
        ifd=2
       endif
       if(dtrend=='d' .and. ifd>0) then
        write(io,'('' You can not detrend the data using both an OLS fit
     * & differencing''/)')
        write(6,'('' You can not detrend the data using both an OLS fit
     *& differencing''/)')
        stop
       endif
      endif
      ts(0)=0.
      ts(1)=0.
      ts(2)=0.
      if(mod=='t') then
c  Linear trend
       par='rcept'
c  Intercept
       ts(0)=rnumread(par,name,delim,ia,ier,io)
c  Slope of the linear part of the trend
       par='lope'
       ts(1)=rnumread(par,name,delim,ia,ier,io)
       par='adratic'
       ts(2)=rnumread(par,name,delim,ia,ier,io)
      endif
      if(mod=='s') then
c  Stochastic trend drift
         par='rift'
         ts(1)=rnumread(par,name,delim,ia,ier,io)
      endif
      if(mod=='t') then
c  Parameters for the curvilinear trend
       par='rcept'
       ts(0)=rnumread(par,name,delim,ia,ier,io)
       if(ier>0) then
        write(io,'(/'' You must set the intercept for the trend'')')
        write(6,'(/'' You must set the intercept for the trend'')')
        stop
       endif
       par='lope'
       ts(1)=rnumread(par,name,delim,ia,ier,io)
       if(ier>0) then
        write(io,'(/'' You must set the linear slope of the trend'')')
        write(6,'(/'' You must set the linear slope of the trend'')')
        stop
       endif
       par='adratic'
       ts(2)=rnumread(par,name,delim,ia,ier,io)
       if(ier>0) then
        write(io,'(/7x,''You must set the quadratic slope of the trend''
     *)')
        write(6,'(/'' You must set the quadratic slope of the trend'')')
        stop
       endif
      endif
c  Write top of ar.out
      if(mod=='t') then
       if(p>0) then
        write(io,'(7x,''Trend plus AR('',i2,'') Model''/)')p
        write(io,'(5x,''Intercept = '',f7.2,2x,''Linear Trend Coefficien
     *t = '',f7.3,/7x,''Quadratic Slope Coefficient = '',f7.3/)')ts(0),t
     *s(1),ts(2)
         elseif(q>0) then
          write(io,'(7x,''Trend plus MA('',i2,'') Model''/)')q
          write(io,'(5x,''Intercept = '',f7.2,2x,''Linear Trend Coeffici
     *ent = '',f7.3,/7x,''Quadratic Slope Coefficient = '',f7.3/)')ts(0)
     *,ts(1),ts(2)
       endif
        elseif(mod=='s') then
         if(p>0) then
          write(io,'(7x,''Stochastic Trend AR('',i2'') with a Drift = ''
     *,f7.2/)')p,ts(1)
         endif
          elseif(mod=='n') then
           if(p>0) then
            write(io,'(7x,''AR('',i2'') Model with no Trend''/)')p
             elseif(q>0) then
              write(io,'(7x,''MA('',i2,'') Model with no Trend''/)')q
           endif
            else
             write(io,*) 'Model switch ',mod,' is incorrect'
             write(6,*) 'Model switch ',mod,' is incorrect'
             stop
      endif
      write(io,'(2x,''Lag for AR fit = '',i6,2x,''AR fit p-value thresho
     *ld'',f7.4/)')np,pv
c-------
c  Forecast range
      nfor=50
c  Burnin no.
      nshft=nint(-4./(alog10(max(rmd,0.9)-1.e-5)))
      ne=n+nshft+nfor+p+1
      na=p-1
      nn=max(p,np)
      nf=max(3*nfor-1,na)
c  Allocate space
      allocate(a(0:nn),x(0:ne),s(0:nf),e(0:ne+n),r(0:ne),c(0:nn),v(0:np*
     *np-1),in(0:nn),mag(0:na),freq(0:na),root(0:na),stat=ibad)
      if(ibad/=0) then
       write(io,*)'Unable to allocate work space for arrays'
       stop
      endif
c  Initialize r & x
      do k=0,ne
       r(k)=0.
       x(k)=0.
      enddo
      ct=0.
c  Initialize a
      do k=0,max(p,np)
       a(k)=0.0
      enddo
c-------
c  Computation of AR or MA parameters
      if(p>1.or.q>1) then
c  AR(p) or MA(q) roots
       pi2=2.d0*dcos(-1.d0)
       do k=0,p/2-1
        mag(k)=dble(rmd)*(1.d0-dfloat(k)/(0.75d0*dfloat(p)))
        freq(k)=pi2*(dble(fr)/dfloat(k+1))
        mag(k+p/2)=mag(k)
        freq(k+p/2)=freq(k)
        root(k)=dcmplx(mag(k)*dcos(freq(k)),-mag(k)*dsin(freq(k)))
        root(k+p/2)=dcmplx(mag(k)*dcos(freq(k)),mag(k)*dsin(freq(k)))
       enddo
       if(p/=2*(p/2)) then
c  p odd
        mag(p-1)=mag(p/2-1)/3.d0
        freq(p-1)=1.d27
        root(p-1)=dcmplx(mag(p-1),0.d0)
       endif
       call polycof(root,p,a,io)
      endif
      if(p==1.or.q==1) then
c  AR(1) or MA(1)
       a(0)=-rmd
       mag(0)=0.d0
       freq(0)=0.d0
      endif
c-------
      if(run==1.and.sig==0.0) then
c  Output impulse return to file ar.data
       do t=0,p
        x(t)=1.0
       enddo
       call argen(r,e,x,a,ct,p,n,ne,0,n,sig,io)
       write(io,'(7x,''Root Max Magnitude = '',f9.3)')rmd
       write(2,'(7x,''Root Max Magnitude = '',f9.3)')rmd
       write(io,'(7x,''1st Root Frequency = '',f5.3/)')fr
       write(2,'(7x,''1st Root Frequency = '',f5.3)')fr
       write(io,'(20(:,'' a('',i2,'')= '',f9.5))')(k+1,-a(k),k=0,p-1)
       if(p>1) then
        write(io,'(/7x,''Amplitudes & Periods of the Roots of the AR Pol
     *ynomial''/)')
        write(io,'(70(:,i3,'' - '',f8.3,2x))')(k+1,mag(k),k=0,p-1)
        write(io,'(70(:,6x,f8.2,2x))')(1.d0/freq(k)*pi2,k=0,p-1)
       endif
       write(io,*)
       write(2,'(20(:,'' a('',i2,'')= '',f9.5))')(k+1,-a(k),k=0,p-1)
       do t=0,n-1
        write(2,'(f15.7)')x(t)
       enddo
       return
      endif
c-------
c  Initialize generator
      call random_seed
      init='n'
      call l90rand(ne,e,ch,init,1)
c-------
      write(io,'(2x,''N = '',i9,2x,''Max Root Magnitude = '',f9.3/)')n,r
     *md
      write(io,'(2x,''Input sig = '',f7.3,7x,''Noise Model is '',a25)')
     *sig,model
      write(io,'(/2x,i8,'' Recursions to achieve stationarity of the AR
     *process'')')nshft
      if(dtrend=='d' .and. nr==1) then
       write(io,'(/2x,''Data detrended using a least squares fit''/)')
      endif
      if(ifd==1 .and. nr==1) then
       write(io,'(/7x,''First difference used to detrend the data'')')
        elseif(ifd==2 .and. nr==1) then
        write(io,'(/7x,''Second difference used to detrend the data'')')
      endif
      if(ot(1:1)=='w') then
       write(io,'(2x,''Random walk of forecasts is used'')')
      endif
      if(ot(2:2)=='r') then
c  Output to AR.log
       if(p>0) then
        write(7,'(2x,''N = '',i9,2x,''Max Root Magnitude = '',f9.3,'' AR
     *('',i2,'') model'')')n,rmd,p
        write(7,'(2x,''Innovations sig = '',f7.3,7x,''Innovations '',a25
     *)')sig,model
       endif
       if(q>0) then
        write(7,'(2x,''N = '',i9,2x,''Max Root Magnitude = '',f9.3,'' MA
     *('',i2,'') model'')')n,rmd,q
        write(7,'(2x,''Innovations sig = '',f7.3,7x,''Innovations '',a25
     *)') sig,model
       endif
       if(dtrend=='d'.and.nr==1) then
        write(7,'(2x,''Data detrended using a least squares fit''/)')
       endif
       if(ifd==1.and.nr==1) then
        write(7,'(7x,''First difference used'')')
         elseif(ifd==2.and.nr==2) then
          write(7,'(7x,''Second difference used'')')
       endif
      endif
c-------
      if(mod=='s') then
c  Stochastic trend with drift
       call argen(r,e,x,a,ct,p,n,ne,nshft,nfor,sig,io)
       x(0)=r(0)
       do t=1,n
        x(t)=ts(1)+x(t-1)+r(t)
       enddo
      endif
      if(p>0.and.mod/='s') then
c  AR(p) model
       call argen(x,e,r,a,ct,p,n,ne,nshft,nfor,sig,io)
      endif
      if(q>0) then
c  MA model
       do t=p,n+p+nfor-1
        sm=0.
        do k=1,p
         sm=sm+a(k-1)*e(t-k)
        enddo
        x(t-p)=sig*sm
       enddo
      endif
c-------
      rn=float(n)
      r16=1./6.
      center=0.5-1./(2.*rn)
      if(mod=='t') then
c  Add linear plus quadratic trend
       do t=0,n+nfor-1
        tc=float(t+1)/rn
        quadt=tc*tc-tc+r16
        trend=ts(0)+ts(1)*(tc-center)+ts(2)*quadt
        x(t)=trend+x(t)
       enddo
      endif
      if(ot(2:2)=='d') then
       if(run==1) then
c  Output data to file ar.data
        write(2,'(2x,''Innovations sig = '',f7.3,2x,''Model '',a25,2x,''
     *Max Root Magnitude = '',f9.3)')sig,model,rmd
        nused=n+nfor
        write(2,'(i9,2x,''1.0 (f15.7)'')')nused
        do t=0,nused-1
         if(mod=='s') then
          write(2,'(f15.5,2x,f15.2)')x(t),ts(1)*float(t)
           else
            write(2,'(f15.7)')x(t)
         endif
        enddo
        close(2)
       endif
      endif
c-------
c  Initialize intercept & slope of estimated linear trend
      d(0)=0.
      d(1)=0.
c  Save last value
      save=x(n-1)
c  Check for differencing
      if(ifd>0) then
       s1=0.d0
       do t=0,n-1
        e(t)=x(t+1)-x(t)
        s1=s1+dble(e(t))
       enddo
       d(0)=x(0)
       d(1)=sngl(s1/dfloat(n))
       if(ifd==2) then
c  Second difference => x
        do t=0,n-2
         x(t)=e(t+1)-e(t)
        enddo
        x(n-1)=0.
         elseif(ifd==1) then
c  First difference => x
          do t=0,n-1
           x(t)=e(t)
          enddo
       endif
      endif
c-------
      write(io,'(1x,77(''*''))')
      at=0.0
      if(dtrend=='d') then
c  Detrend with trend in array e
       do t=0,n-1
        tc=float(t+1)/rn
        e(t)=tc-center
        e(n+t)=tc*tc-tc+r16
       enddo
c  Ordinary least squares fit of curvilinear trend
       po=2
       call ols(x,e,c,v,r,r2,at,se,n,po,ier,io)
       if(ier>0) then
        write(io,*) 'Error in ols ',ier
        stop
       endif
c  Write R square
       write(io,'(7x,''Adjusted R Square = '',f7.3,2x,''Standard Error o
     *f Detrended Data = '',f10.3/)') r2,se
       write(io,'(7x,''Estimate of Intercept = '',f9.3/)') at
       d(0)=at
c  Write out estimates
       do k=0,po-1
        seb=sqrt(sngl(v(k*po+k)))
        tst=sngl(c(k))/seb
        write(io,'(2x,''b('',i2,'') = '',f7.3,2x,''b^ = '',f7.3,2x,''Std
     * Error(b)) = '',g9.3,2x, ''t = '',f10.3)')k+1,ts(k+1),c(k),seb,tst
       enddo
       write(io,'(1x,77(''-''))')
c  Slopes
       d(1)=sngl(c(0))
       d(2)=sngl(c(1))
c  Detrended residuals => x
       do t=0,n-1
        x(t)=r(t)
       enddo
      endif
c-------
c  Statistics for x or detrended x
      call stat(n,x,am,sd,sk,c4,c6,xmax,xmin,io)
      if(np>0) then
c  AR(np) fit
       call far(x,c,s,r,in,r2,se,n,np,po,pv,ier,io)
       if(ier==1) then
        write(io,*) 'Error in subroutine far ier = ',ier
        stop
       endif
       if(po==0) then
        write(io,'(/'' No lag is significant after the iterative least s
     *quares AR fit is done'')')
        write(io,'(/7x,''po = '',i1,2x,''Run = '',i3/)')po,nr
        return
       endif
c-------
       npo=po-1
       write(io,'(/17x,''AR('',i3,'') parameters / t values'')')po
       write(io,'(/6x,''Adjusted R Square = '',f5.3,2x,''Std Error of AR
     * Fit = '',g10.3)')r2,se
       write(io,'(7x,''p-Value Threshold for the Iterative AR Prewhiteni
     *ng Method = '',f5.3)')pv
c  Write out AR coefficients
       write(io,'(/7x,''Estimates of the AR Coefficients and their Stati
     *stics'')')
       write(io,'(1x,57(''=''))')
       do k=0,npo
        write(io,'(2x,''a('',i2,'') = '',f12.5,2x,''a^ = '',f12.5,2x,''S
     *td Error(a^) '',f9.4,2x,''t= '',f11.2)')in(k),a(in(k)-1),-c(k),s(k
     *),-sngl(c(k))/s(k)
       enddo
       if(p>1) then
        write(io,'(1x,57(''=''))')
        write(io,'(/7x,''True Amplitudes & Periods of the Roots of the A
     *R Polynomial''/)')
        write(io,'(70(:,i3,'' - '',f8.3,2x))')(k+1,mag(k),k=0,p-1)
        write(io,'(70(:,6x,f8.2,2x))')(pi2/freq(k),k=0,p-1)
       endif
       write(io,'(1x,57(''=''))')
       write(io,'(/7x,''Amplitudes & Periods of the Roots of the AR Poly
     *nomial'')')
       write(io,'(11x,''A real root is given a period = 0.'')')
       call out(c,in,nn+1,po,io)
       if(mod=='s'.and.ifd>0) then
c  Estimate of drift in the stochastic trend AR model from 1st differences
        write(io,'(1x,77(''*''))')
        write(io,'(7x,''Drift Estimate = '',f10.4)')d(1)
       endif
c-------
c  Statistics
       write(io,'(1x,57(''=''))')
       if(dtrend=='n' .and. ifd==0) then
        write(io,'(/7x,''Statistics for the AR series'')')
        write(io,'(1x,57(''*''))')
         elseif(dtrend=='d') then
          write(io,'(/7x,''Statistics for the detrended series''/)')
          write(io,'(1x,57(''*'')/)')
           elseif(ifd>0) then
            write(io,'(/7x,''Statistics for the differenced series'')')
            write(io,'(1x,57(''*'')/)')
       endif
       write(io, '(7x,'' Mean = '',g12.4,2x,''Sig = '',g12.4/)')am,sd
       write(io,'(7x,'' Skewness = '',g12.3,2x,''Kurtosis = '',g14.2/)')
     *sk,c4
       write(io,'(7x,'' Max ='',g17.5,2x,''Min ='',g17.5)') xmax,xmin
       write(io,'(/1x,57(''=''))')
c  Statistics of resids r <= far
       ip=po+1
       call stat(n-po,r(ip),rm,rs,rk,c4r,c6,rmax,rmin,io)
       exp=0.4
c  Standardize r
       do k=ip,n-po
        r(k)=(r(k)-rm)/rs
       enddo
       call port(r(ip),n-po,exp,ct,io)
       write(io,'(/7x,''Whiteness test statistic p-value = '',f7.3)')ct
       write(io,'(/7x,''Statistics for the Residuals'')')
       write(io,'(1x,57(''*'')/)')
       write(io, '(7x,'' Mean = '',g12.4,2x,''Sig = '',g12.4/)') rm,rs
       write(io,'(7x,'' Skewness = '',g12.3,2x,''Kurtosis = '',f12.2/)')
     *rk,c4r
       write(io,'(7x,'' Max ='',g17.5,2x,''Min ='',g17.5/)') rmax,rmin
c--------
       if(ot(2:2)=='r') then
c  Output to AR.log
        do k=0,npo
         write(7,'(2x,''a('',i2,'') = '',f7.3,2x,''a^ = '',f7.3,2x,''Std
     * Error(a^) '',f9.4,2x,''t= '',f11.2)')in(k),a(in(k+1)),-c(k),s(k),
     *-sngl(c(k))/s(k)
        enddo
        call out(c,in,np,po,7)
        write(7,'(1x,57(''*'')/)')
       endif
      endif
c  End of AR fit
c-------
      if(ot(1:1)/='n'.and.run==1) then
       write(4,'(10x,''x(t)'',5x,''Forecast - True a'',3x,''Forecast - a
     *^'')')
       sg=0.
c  Shift back last block of x's
       do t=0,p-1
        r(t)=x(n-p+t)
       enddo
       if(ot(1:1)=='f') then
c  Forecast AR using true parameters
        call argen(s,e,r,a,am,2,nfor,ne,nshft,nfor,sg,io)
c  Forecast AR using estimated parameters
        do t=0,npo
         a(t+1)=-sngl(c(t))
         r(t)=x(n-po+t)
        enddo
        call argen(s(nfor),e,r,a,am,po,nfor,ne,nshft,nfor,sg,io)
       endif
c  Output true x's for t=n,n+nfor-1 & forecasts
       trend=0.
       do t=0,nfor-1
        if(dtrend=='d' .and. ot(1:1)=='f') then
c  Forecast with trend & AR
         tc=float(n+t+1)/rn
         quadt=tc*tc-tc+r16
         trend=d(0)+d(1)*(tc-center)+d(2)*quadt
         write(4,'(3(f15.4,1x))')x(n+t),s(t)+trend,s(nfor+t)+trend
        endif
        if(ifd==1 .and. ot(1:1)=='f') then
c  Forecast with line from 1st difference & AR
         trend=d(0)+d(1)*float(n+t+1)
         write(4,'(3(f15.4,1x))')x(n+t),s(t)+trend,s(nfor+t)+trend
        endif
c  Forecast AR
        if(dtrend=='n'.and.ot(1:1)=='f') then
         write(4,'(3(f15.4,1x))')x(n+t),s(t),s(nfor+t)
        endif
c  Forecast with mean drift & AR
        if(ifd==1.and.ot(1:1)=='w') then
         write(4,'(2(f15.4,1x))')x(n+t),save+d(1)*float(t+1)
        endif
       enddo
       call flush(4)
       close(4)
      endif
c-------
      deallocate(a,x,s,e,r,c,v,in,mag,freq,root)
      return
      end
c
c============================================
c
      subroutine argen(x,e,r,a,ct,p,n,ne,nshft,nfor,sig,io)
c********************************************
c  If sig>0 n+nfor values AR(p) x(t)+a(1)*x(t-1)+...+a(p)*x(t-p)=ct+e(t) with
c   burnin nshft
c  If sig=0 then x(t+1)=ct-a(1)*x(t)+...-a(p)*x(t-p)
c  Initial values - r(1),...,r(p)
c============================================
      integer p,t
      real x(n+nfor),e(ne),r(ne),a(p)
      intent(in) a,e,ct,p,n,ne,nshft,nfor,sig,io
      intent(out) x
      intent(in out) r
c********************************************
c  Trap p
      if(p>50) then
       write(io,*)'p ',p,' ARGEN > 50'
       return
      endif
      if(sig>0.0) then
c  AR(p)
       nsh=nshft
       nuse=n+nfor
      endif
      if(sig<=0.) then
c  Convolution with no noise
       nsh=p
       nuse=nfor
      endif
c  Start recursions
      do t=1,nuse+nsh
       sm=0.
       do k=1,p
        sm=sm-a(k)*r(t+p-k)
       enddo
       if(sig>0.0) then
        r(t+p)=ct+sm+sig*e(t)
       endif
       if(sig<=0.0) then
        r(t+p)=ct+sm
       endif
      enddo
      do t=1,nuse
       x(t)=r(t+nsh)
      enddo
c-------
      return
      end
c
c===========================================
c
      subroutine ols(y,x,c,a,r,r2,am,se,n,p,ier,io)
c********************************************
c  Ordinary least squares fit of y=cx+e. Version 6-24-00
c  Iput: x - array of nxp data matrix in column major order
c    y - dependent variable observations
c  Output: c - estimates of c coefficients, am - estimate of intercept
c    r(1),...,r(n) - residuals, a - covariance matrix of estimates
c    r2 - Adjusted R square, se - stnd. deviation of residuals
c
c   ier = 1: A is not positive definite, ier = 2: p < 1,
c   ier = 3: p > n/2, ier = 4: work array can not be allocated
c============================================
      allocatable::v,b,w,d,m
      integer p
      real x(n*p),y(n),r(n),m(:),v(:),w(:)
      real*8 c(p),a(p*p),b(:),d(:)
      intent(in) x,y,n,p,io
      intent(out) c,a,r,r2,am,se,ier
c********************************************
c  Trap p
      if(p<1) then
       write(io,*) 'p in OLS(p) ',p,' < 1'
       ier=2
       return
        elseif(p>n/2) then
         write(io,*) 'p in OLS(p) ',p,' > n/2 ',n/2,' n = ',n
         ier=3
         return
      endif
c-------
      na=p*p-1
      np=p-1
c  Allocate space
      allocate(v(0:na),b(0:np),w(0:np),d(0:np),m(0:np),stat=ibad)
      if(ibad/=0) then
       write(io,*) 'Unable to allocate work space for work array in ols'
       ier=4
       return
      endif
      call ls(y,x,c,a,r,v,b,w,d,m,n,p,r2,am,se,ier,io)
c-------
      deallocate(v,b,w,d,m)
      return
      end
c
c============================================
c
      subroutine ls(y,x,c,a,r,v,b,w,d,m,n,p,r2,am,se,ier,io)
c********************************************
c  Least squares (OLS) fit of y=bx+e
c============================================
      integer p,t
      real y(n),x(n*p),m(p),r(n),v(p*p),w(p)
      real*8 c(p),a(p*p),b(p),d(p),dy,sy,ssy,sse,sm
      intent(in) x,y,n,p,io
      intent(out) c,a,r,r2,am,se,ier
c********************************************
      call cva(y,x,a,v,b,w,m,sy,ssy,n,p,io)
c  Cholesky decomposition of a
      call dcholdc(a,p,d,ier,io)
      if(ier/=0) then
       write(io,*) 'DCHOLDC failed, A is not positive definite'
       ier=1
       return
      endif
      call dcholsl(a,p,d,b,c)
c  Inverse of A
      call dinverse(a,p,d)
c-------
c  Residuals
      sse=0.d0
      do t=1,n
       sm=0.d0
       do k=1,p
        kt=(k-1)*n+t
        sm=sm+c(k)*dble(x(kt)-m(k))
       enddo
       dy=dble(y(t))-sy-sm
       sse=sse+dy*dy
       r(t)=sngl(dy)
      enddo
c  Standard error of residuals
      se=sqrt(sngl(sse)/float(n-p))
c  Intercept
      cn=0.
      do k=1,p
       cn=cn+sngl(c(k))*m(k)
      enddo
      am=sngl(sy)-cn
c  R squared
      r2=1.-sngl(sse/ssy)
      r2=r2-float(p)/float(n-p-1)*(1-r2)
c  Covariance matrix of estimates in a
      do j=1,p
       do k=j,p
        a((k-1)*p+j)=se*a((k-1)*p+j)
        a((j-1)*p+k)=a((k-1)*p+j)
       enddo
      enddo
c-------
      return
      end
c
c===========================================
c
      subroutine cva(y,x,a,v,b,w,m,sy,ssy,n,p,io)
c********************************************
c  Means and covariances of x & y
c============================================
      integer p,t
      real y(n),x(n*p),m(p),v(p*p),w(p)
      real*8 a(p*p),b(p),sy,ssy,y1,y2,sx,sb,rn
      intent(in) x,y,n,p,io
      intent(out) a,b,m,v
c********************************************
      rn=dfloat(n)
c  Compute means
      y1=0.d0
      y2=0.d0
      do k=1,p
       sx=0.d0
       do t=1,n
        dy=dble(y(t))
        if(k==1) then
         y1=y1+dy
         y2=y2+dy*dy
        endif
        kt=(k-1)*n+t
        sx=sx+dble(x(kt))
       enddo
c  Means of the x's
       m(k)=sngl(sx/rn)
      enddo
c-------
c  Mean of y
      sy=y1/rn
c  Sum of squares of y-mean(y)
      ssy=y2-rn*sy*sy
c  Covariances in upper triangle of matrix a
      do j=1,p
       do k=j,p
        sx=0.d0
        sb=0.d0
        do t=1,n
         jt=(j-1)*n+t
         kt=(k-1)*n+t
         sx=sx+dble((x(jt)-m(j))*(x(kt)-m(k)))
         if(k==p) then
          sb=sb+(dble(y(t))-sy)*dble(x(jt)-m(j))
         endif
        enddo
        a((k-1)*p+j)=sx
        v((k-1)*p+j)=sngl(sx)
       enddo
       b(j)=sb
       w(j)=sngl(sb)
      enddo
c-------
      return
      end
c
c===========================================
c
      subroutine out(c,in,maxp,po,io)
c*******************************************
c  Computes the roots of the polynomial whose coefficients are in c
c===========================================
      allocatable::a,rt
      integer po,in(maxp),ep
      real*8 c(maxp),drer,dimr,per,pi
      complex a(:),rt(:)
      logical*1 plsh
      intent(in) c,maxp,po,io
c********************
      allocate(a(0:maxp),rt(0:maxp-1),stat=ibad)
      if(ibad/=0) then
       write(io,*) 'Unable to allocate work space for arrays in out'
       stop
      endif
      ep=in(po)
      do k=0,maxp
       a(k)=cmplx(0.,0.)
      enddo
c-------
      pi=dcos(-1.d0)
      plsh=.true.
c  Put coefficients of AR into a
      do k=1,ep
       a(ep-in(k))=cmplx(-1.*sngl(c(k)),0.)
      enddo
      ier=0
      a(ep)=cmplx(1.,0.)
      call zroots(a,ep,rt,plsh,ier)
      if(ier==1) then
       write(io,'(/7x,''No convergence in AR roots calculation''/)')
       return
      endif
c-------
      do k=0,ep-1
       amp=cabs(rt(k))
       if(abs(aimag(rt(k)))<1.e-15) then
        per=0.d0
         else
          drer=dble(real(rt(k)))
          dimr=dble(aimag(rt(k)))
          per=(2.d0*pi)/datan2(dimr,drer)
          if(per>1.d+4) then
           per=1.d4
          endif
       endif
       rt(k)=cmplx(amp,sngl(abs(per)))
      enddo
c-------
      write(io,'(1x,57(''=''))')
      write(io,'(70(:,i3,'' - '',f8.3,2x))') (k+1,real(rt(k)),k=0,ep-1)
      write(io,'(70(:,6x,f8.2,2x))') (aimag(rt(k)),k=0,ep-1)
c-------
      return
      end