FRACTINT.RPF

FRACTINT.RPF is an example of estimation of an ARFIMA Model using frequency domain methods.

An ARFIMA model takes the form

\begin{equation} \left( {1 - L} \right)^d A\left( L \right)y_t = B\left( L \right)u_t \end{equation}

where $A$ and $B$ are standard autoregressive and moving average polynomials. Solving this out gives

\begin{equation} y_t = \left( {{{B\left( L \right)} \mathord{\left/ {\vphantom {{B\left( L \right)} {A\left( L \right)}}} \right. } {A\left( L \right)}}} \right)\left( {1 - L} \right)^{ - d} u_t \end{equation}

The square of the transfer function of the lag polynomial on the right is the spectral density function that we want (except for scale by the variance of $u$). The easiest way to compute this is to use a complex-valued FRML to give us the transfer function itself. This makes it simpler to change from one model to another, since the model is entirely described by the transfer function.

If you use this method of estimation, there are two things you must be sure to do:

1.Run the MAXIMIZE instruction across the full range of frequency domain ordinates. This is, almost certainly, longer than your actual data set.

2.Include the extra factor of actual observations divided by the number of frequency ordinates. Because of the procedure’s reliance on frequency domain filtering, it only works acceptably if the data are padded. This factor corrects for the effect of padding.

As described above, this starts out with a complex-valued FRML, which computes the transfer function. The chosen model is an ARFIMA(1,d,1) which has four free parameters: the autoregressive parameter (A), the moving average parameter (B), the fractional difference parameter (D) and the innovation variance (IVAR). In the FRML, 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.

declare frml[complex] transfer

nonlin a b d ivar

frml transfer = (1+b*%zlag(t,1))/$

((1-a*%zlag(t,1))*%conjg((1-%zlag(t,1))^D))

This next runs an AR(1) solely to get a reasonable guess for the variance. None of the other parameters are scale-sensitive, so they’re given basic guesses.

linreg(noprint) usasur

# constant usasur{1}

compute a=.9, b=0.0, d=.5, ivar=%seesq

Following the standard practice for frequency-domain filter, this uses double the recommended number of ordinates,

compute nobs = %allocend()

compute nords = 2*%freqsize(nobs)

freq 3 nords

For the so-called Whittle likelihood, the (unsmoothed) periodogram provides sufficient statistics for the data. We need a real-valued copy of that since MAXIMIZE eventually needs a real-valued log likelihood.

rtoc

# unemp0

# 1

fft 1

cmult(scale=1.0/nobs) 1 1 / 2

set periodogram 1 nords = %real(%z(t,2))

The log likelihood for the sample is computed by summing log likelihood elements for each ordinate. Given a particular set of parameters, the theoretical spectral density is the variance times the square of transfer function. The observed “data” value is the periodogram. The %LOGDENSITYCV function uses these to compute the log likelihood of the ordinate. Because there are NORDS ordinates, but only NOBS observations, the number of observations parameter in the %LOGDENSITYCV is set to the fraction that corrects for that. MAXIMIZE then can estimate the model as with any other log likelihood, except that it needs to be summed across the number of ordinates, rather than the number of observations.

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 Fourier Transform of the series by the transfer function gives the Fourier Transform of the residuals. Inverse transform and send the relevant part back to the time domain. The first few residuals in this as they are likely to be very large with a correspondingly large variance, so they are skipped from the residual autocorrelation analysis:

cset 3 = %z(t,1)/transfer(t)

ift 3

ctor 1 nobs

# 3

# resids

corr(qstats,span=24) resids 3 *

This computes an graphs the corresponding Moving Average Representation (MAR) for the model.

set mar 1 100 = %real(%z(t,3))
graph(style=bargraph,number=0,header="Moving Average Representation")
# mar 1 100

Full Program

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

Output

MAXIMIZE - Estimation by BFGS

Convergence in 9 Iterations. Final criterion was 0.0000021 <= 0.0000100

Quarterly Data From 1960:01 To 2087:04

Usable Observations 511

Skipped/Missing (from 512) 1

Function Value 286.6385

Variable Coeff Std Error T-Stat Signif

************************************************************************************

1. A 0.7668430118 0.1134879584 6.75704 0.00000000

2. B 0.1216759954 0.1201787621 1.01246 0.31131893

3. D 0.6276904559 0.2049611620 3.06248 0.00219508

4. IVAR 0.0029101357 0.0002997199 9.70952 0.00000000

Correlations of Series RESIDS

Quarterly Data From 1960:03 To 2007:02

Autocorrelations

1 2 3 4 5 6 7 8 9 10

0.14489 0.06591 0.07627 -0.12325 -0.10977 0.06652 0.08206 -0.19079 -0.01747 0.03310

11 12 13 14 15 16 17 18 19 20

-0.08988 -0.05656 -0.00342 -0.07343 -0.01969 -0.01437 0.00394 -0.02350 0.04854 -0.01726

21 22 23 24 25 26 27

0.05314 0.01970 -0.07772 -0.00827 0.00047 -0.00879 -0.08541

Ljung-Box Q-Statistics

Lags Statistic Signif Lvl

24 27.172 0.296487

Graph