*
* FRACTINT.RPF
* RATS Version 8, User's Guide, Example 14.3
*
open data oecdsample.rat
calendar(q) 1960
data(format=rats) 1960:1 2007:2 usaunempsurqs
*
* The data series is the US unemployment rate. The mean is removed after
* transformation to square roots.
*
set usasur = sqrt(usaunempsurqs)
diff(center) usasur / unemp0

declare frml[complex] transfer
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*%zlag(t,1))/$
   ((1-a*%zlag(t,1))*%conjg((1-%zlag(t,1))^D))
*
* Run an AR(1) to get a guess for the innovation variance.
*
linreg(noprint) usasur
# constant usasur{1}
*
compute a=.9,b=0.0,d=.5,ivar=%seesq
*
* Use double the recommended number of ordinates. Send the data to the
* frequency domain and compute the periodogram.
*
compute nobs  = %allocend()
compute nords = 2*%freqsize(nobs)
freq 3 nords
rtoc
# unemp0
# 1
fft 1
cmult(scale=1.0/nobs) 1 1 / 2
set periodogram 1 nords = %real(%z(t,2))
frml likely    = trlambda=transfer, $
  %logdensitycv(ivar*%cabs(trlambda)**2,periodogram,(float(nobs)/nords))
maximize(pmethod=simplex,piters=5,method=bfgs) 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 1 nobs
# 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 3 *
*
* Compute the (approximate) moving average representation
*
* The transfer function doesn't exist at zero frequency for d>0. (The
* 0^0 calculation in computing "transfer" produces a 0). Since the MAR
* is known to have a unit lead coefficient, we know what we need to add
* to the inverse transform to achieve that. We just add that same
* adjustment term to all the lags to produce (approximately) the correct
* result.
*
cset 3 1 nords = transfer
ift 3
compute adjust=1-%z(1,3)
cset 3 1 nords = %z(t,3)+adjust
*
set mar 1 100 = %real(%z(t,3))
graph(style=bargraph,number=0,header="Moving Average Representation")
# mar 1 100