* ARFSIM.SRC
*
* Frequency domain procedure for generating a true ARFIMA(0,d,0)
*
* Programmed by Rob Schoen (rjschoen@earthlink.net or rjs2@stern.nyu.edu)
* December 19, 1996
*
* Davies and Harte Method (see Beran p 216-217)
* Simulates ARFIMA(0,d,0) model with unit innovation variance.
* Gam is the true autocovariance sequence
*
* (If you've been using FIF.SRC (fractional integration filter)
* to simulate long memory, you will find this procedure to be orders
* of magnitude faster. It also produces authentic ARFIMA's, while
* FIF is just an approximation.)
*
* Syntax:
* @ARFSIM d n series
* d is the fractional differencing parameter (works for -1 <= d < 1.5)
* n is the number of observations to be generated
* series is an empty series into which the ARFIMA realization is placed.
*
*
* Critical Info:
*
* There must be room allocated for 2*n observations in all real series.
*
* A frequency instruction must be run prior to 'sourcing' this procedure
* (see ARFSIM.PRG).
*
* This procedure requires the XGAMMA procedure, which must be 'sourced'
* before this procedure is 'sourced'(see ARFSIM.PRG). Rats doesn't have
* a gamma function; you have to use exp(%lngamma(x)), but this doesn't
* work for x < 0. XGAMMA will calculate gamma for -1 < x < 0, which is
* necessary for this procedure.


procedure arfsim d n x
type real *d
type integer *n
type series *x

declare cseries cb z

local integer Noninvert Nonsta r h
local real gx g1 g2 g3 g4 g5 gz
local series Gam U V gg xx fu fv

clear x

compute Noninvert = 0
compute Nonsta = 0


* This section checks if the ARFIMA requested is stationary and invertible.
* If it's non-stationary, d-1 will be simulated and then integrated.
* If it's non-invertible, d+1 will be simulated and then differenced.

if (d>=-1).and.(d<-.5)
   {
   compute Noninvert = 1
   compute d = d + 1
   compute n = n + 1
   }
if (d>=.5).and.(d<1.5)
   {
   compute Nonsta = 1
   compute d = d - 1
   }
if d==0
   {
   set x 1 n = %ran(1)
   }

if d<>0
   {

   * This section calculates the true autocovariance sequence

   clear Gam U V gg

   compute r = 0
   @xgamma (r+d) g1
   @xgamma (1.0-2*d) g2
   @xgamma (r-d+1.0) g3
   @xgamma (d) g4
   @xgamma (1.0-d) g5

   compute Gam(1) = g1*g2/(g3*g4*g5)

   @xgamma d g4

   do r = 1,n-1

	  compute g1 = (%lngamma(r+d))
	  compute g2 = (%lngamma(1.0-2*d))
	  compute g3 = (%lngamma(r-d+1))
	  compute g5 = (%lngamma(1-d))

	  compute Gam(r+1) = g1+g2-(g3+g5)
	  compute Gam(r+1) = exp(Gam(r+1))/g4

   end do r

   do r = 1,n,1
      compute gg(r) = Gam(r)
   end do r
   compute h = n
   do r = n-1,2,-1
      compute h = h + 1
      compute gg(h) = Gam(r)
   end do r

   rtoc
   # gg
   # cb
   fft cb 1 h
   conjugate cb 1 h

   * This section generates two series of standard normals

   set u 1 n = %ran(1)
   compute u(1) = sqrt(2.0)*u(1)
   compute u(n) = sqrt(2.0)*u(n)
   set v 1 n = %ran(1)
   compute v(1) = 0.0
   compute v(n) = 0.0

   clear fu fv xx
   do r = 1,n
      compute fu(r) = u(r)
      compute fv(r) = v(r)
   end do r
   compute h = n
   do r = n-1,2,-1
      compute h = h + 1
      compute fu(h) = u(r)
      compute fv(h) = -v(r)
   end do r

   * This section uses FFT to generate the ARFIMA series.

   cset z 1 h = %cmplx(fu,fv)

   csqrt cb 1 h
   cset z 1 h = (cb*z)
   fft z 1 h
   conjugate z 1 h
   cset z 1 h = 1/(2*sqrt(n-1))*z
   ctor
   # z
   # xx
   set x 1 n = xx
   }

* If the request was for a noninvertible ARFIMA, it is differenced.
* If the request was for a nonstationary ARFIMA, it is integrated.

if noninvert==1
   {
   diff x
   compute n = n - 1
   compute d = d - 1
   }
if nonsta==1
   {
   accumulate x
   compute d = d + 1
   }

END

