* * FRACTINT.RPF * Estimation of an ARFIMA Model * * RATS User's Guide, example from Section 14.10 * 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 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