Preamble

library("forecast")
library("fracdiff") # for fracdiff, fdGPH
library("tseries")  # for garch

Data Sets

# US CPI, Jan 1947 - Jan 2017, Monthly
cpi.data <- read.csv("http://ptrckprry.com/course/forecasting/data/cpi.csv")
cpi <- cpi.data$cpi
log.cpi <- log(cpi)
# US Unemployment Rate, Jan 1948 - Jan 2017, Montly, Seasonally adjusted
unemp.data <- read.csv("http://ptrckprry.com/course/forecasting/data/unemployment.csv")
unemp <- unemp.data$unemployment
log.unemp <- log(unemp)
# Russell 2000 Index, Sep 10, 1987 - Feb 10, 2017, Daily
russell.data <- read.csv("http://ptrckprry.com/course/forecasting/data/russell.csv")
russell <- russell.data$russell
log.russell <- log(russell)
# VIX Volatility Index, Jan 2, 2004 - Apr 20, 2017, Daily
vix.data <- read.csv("http://ptrckprry.com/course/forecasting/data/vix.csv")
vix <- vix.data$vix.adjclose
log.vix <- log(vix)
logrange.vix <- log(vix.data$vix.hi) - log(vix.data$vix.lo)
# SNP500 Daily Realized Volatility,  Feb 1, 1983 - June 29, 2000
snpvol.data <- read.csv("http://ptrckprry.com/course/forecasting/data/SNPVol.csv")
snpvol <- snpvol.data$SNPVol
log.snpvol <- log(snpvol)
make.plots <- function(x, main="", lag.max=100)
{
    ylab <- deparse(substitute(x))

    fit <- fdGPH(x)

    par(mfrow=c(1,2))
    plot(x, type="l", xlab="time", ylab=ylab, col=2, main=main)
    mtext(paste0("d = ", round(fit$d, 2), ", SE = ", round(fit$sd.reg, 2)))
    Acf(x, lag.max=lag.max, main=paste0("Series: ", ylab))
}

make.plots(log.russell, "Log(Russell 2000 Index), Sep 10, 1987 - Feb 10, 2017")

make.plots(log.cpi, "Log(US CPI), Jan 1947 - Jan 2017")

make.plots(log.vix, "Log(VIX Spot Volatility Index), Jan 2, 2004 - Apr 20, 2017")

make.plots(log.unemp, "Log(US Unemployment), Jan 1948 - Jan 2017")

make.plots(log.snpvol, "Log(Daily Realized S&P500 Vol), Feb 1, 1983 - June 29, 2000")

make.plots(logrange.vix, "VIX Log Range, Jan 2, 2004 - Apr 20, 2017")

Dow

# DJIA, Dec 12, 1914 - Feb 6, 2017
data.dow <- read.csv("http://ptrckprry.com/course/forecasting/data/dow.csv")
dow <- data.dow$dow

Log(Dow)

log.dow <- log(dow)
plot(log.dow, type="l", col=2)

Dow Returns (Diff Log)

ret.dow <- diff(log.dow)
plot(ret.dow, type="l", col=2)

Long Memory in Dow Volatility (Squared Returns)

plot(ret.dow^2, type="l", col=2)

Acf(ret.dow^2, lag.max=200)

Pacf(ret.dow^2, lag.max=200)

GARCH(1,1) Fit

garch(ret.dow, c(1,1), trace=FALSE)

Call:
garch(x = ret.dow, order = c(1, 1), trace = FALSE)

Coefficient(s):
       a0         a1         b1  
1.019e-06  7.854e-02  9.133e-01

Dow Returns: d = 0 or d = 1?

Returns

plot(ret.dow, type="l", col=2)

par(mfrow=c(1,2))
Acf(ret.dow, lag.max=200)
Pacf(ret.dow, lag.max=200)

Diff. Returns

plot(diff(ret.dow), type="l", col=2)

par(mfrow=c(1,2))
Acf(diff(ret.dow), lag.max=200)
Pacf(diff(ret.dow), lag.max=200)

Methods to Estimate d

Method 1

fit1 <- fracdiff(ret.dow)
summary(fit1) # SE is probably bogus here

Call:
  fracdiff(x = ret.dow) 

Coefficients:
   Estimate Std. Error z value Pr(>|z|)    
d 0.0107638  0.0008766   12.28   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
sigma[eps] = 0.01098304 
[d.tol = 0.0001221, M = 100, h = 0.0008918]
Log likelihood: 8.464e+04 ==> AIC = -169267.3 [2 deg.freedom]

Method 2

fit2 <- fdGPH(ret.dow) # more reliable
cat("d: ", fit2$d, "\n")
d:  0.04301893 
cat("SE: ", fit2$sd.reg, "\n")
SE:  0.05271393

Method 3

fit3 <- arfima(ret.dow) # sometimes works, sometimes crashes R
summary(fit3)

Call:
  arfima(y = ret.dow) 

*** Warning during (fdcov) fit: unable to compute correlation matrix; maybe change 'h'

Coefficients:
       Estimate
d         0.191
ar.ar1    0.927
ar.ar2    0.169
ar.ar3   -0.142
ma.ma1    1.091
ma.ma2    0.132
ma.ma3   -0.242
sigma[eps] = 0.01097119 
[d.tol = 0.0001221, M = 100, h = 0.0008921]
Log likelihood: 8.467e+04 ==> AIC = -169314.1 [8 deg.freedom]

Simulating ARFIMA

x <- fracdiff.sim(n = 5000, d = 0.0164)$series

TS Plot

plot(x, type="l", col=3)

ACF, PACF

par(mfrow=c(1,2))
Acf(x, lag.max=200)
Pacf(x, lag.max=200)