*
* Example 10.5.1 on pp 363-365
* (Note-this takes several minutes to run)
*
open data nile.tsm
calendar 622
data(format=free,org=columns) 622:01 871:01 nile
*
diff(center) nile / cnile
@bjautofit(pmax=10,qmax=10) cnile
*
boxjenk(maxl,ar=5,ma=3,demean) nile
*
declare frml[complex] transfer
dec vect a(2) b(2)
nonlin a b d ivar
*
* Construct the transfer function of the moving average polynomial. This example
* includes one MA lag and one AR lag in addition to the fractional difference.
* The %conjg function on the final term is necessary in the next function since
* the * operator conjugates the right operand, and we want a straight multiply.
*
frml transfer = (1+b(1)*%zlag(t,1)+b(2)*%zlag(t,2))/$
   ((1-a(1)*%zlag(t,1)-a(2)*%zlag(t,2))*%conjg((1-%zlag(t,1))**D))
*
* Run an AR(1) to get a guess for the innovation variance.
*
linreg(noprint) nile
# constant nile{1}
*
compute a=||.3,.3||,b=||0.,.5||,d=.5,ivar=%seesq
*
* Use double the recommended number of ordinates. Send the data to the frequency
* domain and compute the periodogram.
*
nonlin a b d=.4 ivar
compute nords = %freqsize(871:1)*2
compute nords = 871:1
compute nobs  = 871:1
freq 3 nords
rtoc
# cnile
# 1
fft 1
cmult(scale=1.0/nobs) 1 1 / 2
frml periodogram = %real(%z(t,2))
frml likely    = trlambda=transfer, (float(nobs)/nords)*$
  (-.5*log(2*%pi)-.5*log(ivar)-log(%cabs(trlambda))-$
     .5/ivar*periodogram/%cabs(trlambda)**2)
maximize(iters=400,method=bfgs,pmethod=simplex,piters=5) likely 1 nords
*
* Dividing the FT of the series by the MAR transfer function gives the FT of the
* residuals. Inverse transform and send the relevant part back to the time domain.
*
cset 3 = %z(t,1)/transfer(t)
ift 3
ctor 622:1 871:1
# 3
# resids
*
* Check out the correlation of residuals. We skip the first few residuals in this
* as they are likely to be very large with a correspondingly large variance.
*
corr(qstats,span=24) resids