library("forecast")
library("fracdiff") # for fracdiff, fdGPH
Raw Data
# US CPI, Jan 1947 - Jan 2017, Monthly
data <- read.csv("http://ptrckprry.com/course/forecasting/data/cpi.csv")
head(data)
date cpi
1 1947-01-01 21.48
2 1947-02-01 21.62
3 1947-03-01 22.00
4 1947-04-01 22.00
5 1947-05-01 21.95
6 1947-06-01 22.08
tail(data)
date cpi
836 2016-08-01 240.389
837 2016-09-01 241.006
838 2016-10-01 241.694
839 2016-11-01 242.199
840 2016-12-01 242.821
841 2017-01-01 244.158
Extract Variables
date <- as.Date(data$date)
cpi <- data$cpi
log.cpi <- log(cpi)
inflation <- c(NA, diff(log.cpi))
Estimate
xbar <- mean(inflation, na.rm=TRUE)
print(xbar)
[1] 0.002893682
Standard Error (?)
n <- sum(!is.na(inflation))
s <- sd(inflation, na.rm=TRUE)
se <- s / sqrt(n)
print(se)
[1] 0.0001194141
95% Confidence Interval (?)
xbar + c(-1,1) * 1.96 * se
[1] 0.002659631 0.003127734
plot(date, inflation, type="l", col=2)
par(mfrow=c(1,2))
Acf(inflation, lag.max=200)
Pacf(inflation, lag.max=200)
par(mfrow=c(1,2))
Acf(diff(inflation), lag.max=200)
Pacf(diff(inflation), lag.max=200)
Method 1
## fit1 <- fracdiff(inflation) # gives an error
fit1 <- fracdiff(inflation[-1]) # remove missing value
summary(fit1)
Call:
fracdiff(x = inflation[-1])
Coefficients:
Estimate Std. Error z value Pr(>|z|)
d 0.3868483 0.0000397 9743 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
sigma[eps] = 0.002723407
[d.tol = 0.0001221, M = 100, h = 3.972e-05]
Log likelihood: 3768 ==> AIC = -7532.733 [2 deg.freedom]
Method 2
fit2 <- fdGPH(inflation[-1])
cat("d: ", fit2$d, "\n")
d: 0.388747
cat("SE: ", fit2$sd.reg, "\n")
SE: 0.118978
Method 3
fit3 <- arfima(inflation)
summary(fit3)
Call:
arfima(y = inflation)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
d 4.583e-05 2.999e-02 0.002 0.999
ar.ar1 9.633e-01 2.868e-03 335.815 <2e-16 ***
ma.ma1 5.450e-01 2.649e-03 205.745 <2e-16 ***
ma.ma2 2.209e-01 2.815e-03 78.476 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
sigma[eps] = 0.002739641
[d.tol = 0.0001221, M = 100, h = 3.967e-05]
Log likelihood: 3764 ==> AIC = -7518.128 [5 deg.freedom]
x <- fracdiff.sim(n = 5000, d = 0.39)$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)
Mean from Simulated Series
(mean.x <- mean(x))
[1] -0.02888058
Naive SE from Simulated Series
sd.x <- sd(x)
n.x <- length(x)
(se.x <- sd.x / sqrt(n.x))
[1] 0.01833684
Naive CI from Simulated Series
mean.x + c(-1,1) * 1.96 * se.x
[1] -0.064820778 0.007059624
Invalid due to Long Memory in Series!
1/sqrt(n) # With indep. observations, SE is of this order
[1] 0.03450328
n^(2 * 0.39 - 1) # With long-memory (d = 0.39), this order instead
[1] 0.2273309