library(ggfortify)
library(forecast)
path_fig = "~/Dropbox/Evry/TimeSeries/Cours/Figures"

An example of ARIMA model

We consider the analysis of quarterly U.S. GNP from 1947(1) to 2002(3), n = 223 observations. The data are real U.S. gross national product in billions of chained 1996 dollars and have been seasonally adjusted. The data were obtained from the Federal Reserve Bank of St. Louis (http://research.stlouisfed.org/). This example is from Peter J Brockwell and Richard A Davis, Time series : theory and methods, Springer Science & Business Media, 2013.

Here are the steps to build an ARIMA model for the gnp series of astsa library.

library(astsa)

1. Time plot of gnp

autoplot(gnp)

There is a clear trend, the variability is perhaps increasing with the mean. We’ll the differentiation on log(gnp) and gnp

ggAcf(gnp)

The ACFs is, as expected, (too) slowly decreasing.

2. Transformation of the data

autoplot(log(gnp))

ggAcf(log(gnp))

3. Differentiation

3.1 For gnp

autoplot(diff(gnp))

ggAcf(diff(gnp))

This is still not stationary: neither at order 1 (mean), nor at order 2 (variance). Let’s try higher orders.

autoplot(diff(gnp,2))

ggAcf(diff(gnp,2))

3.2 For log(gnp)

autoplot(diff(log(gnp)))

ggAcf(diff(log(gnp)))

This looks reasonability stationnary at orders 1 and 2. We now construct an ARMA model for diff(log(gnp))

4. ACFs and PACFs for diff(log(gnp))

ggAcf(diff(log(gnp)))

ggPacf(diff(log(gnp)))

The ACF calls for an MA(2), while PACF calls for an AR(1). We’ll then estimate in ARIMA(1,1,0) and ARIMA(0,1,2) models for log(gnp).

5. Model selection

Caution: the default for the option include.drift in the Arima function of the forecast package is FALSE, but in our case the diff(log(gnp)) is clearly not centered

autoplot(diff(log(gnp)))

arima_110 = Arima(log(gnp),order = c(1,1,0),include.constant=TRUE)
arima_012 = Arima(log(gnp),order = c(0,1,2),include.constant=TRUE)
criteria = matrix(c(arima_110$aicc,arima_012$aicc,arima_110$bic,arima_012$bic),nrow=2,byrow = TRUE)
colnames(criteria) = c("arima_110","arima_012")
rownames(criteria) = c("aicc","bic")
criteria
     arima_110 arima_012
aicc -1431.110 -1431.745
bic  -1421.012 -1418.318

6. Diagnosis on residuals

check = checkresiduals(arima_110 )

    Ljung-Box test

data:  residuals
Q* = 9.8183, df = 6, p-value = 0.1325

Model df: 2.   Total lags used: 8

auto.arima(log(gnp))

An example of SARIMA model

A pure pure seasonal \(ARMA(0,1)_S\) model

set.seed(666)
omega =rnorm(n=1000)
X = ts(omega[13:1000] - 0.5 * omega[1:988], start = c(1900,1), frequency = 12)
autoplot(X)

ggAcf(X)

dev.print(pdf,paste(c(path_fig,"acf_pure_season_MA1.pdf"),collapse="/"))
quartz_off_screen 
                2 
ggPacf(X)

dev.print(pdf,paste(c(path_fig,"pacf_pure_season_MA1.pdf"),collapse="/"))
quartz_off_screen 
                2 
phi = c(rep(0,11),.8)
ACF = ARMAacf(ar=phi, ma=-.5, 50)[-1]     # [-1] removes 0 lag
PACF = ARMAacf(ar=phi, ma=-.5, 50, pacf=TRUE)
par(mfrow=c(1,2), mar=c(2.5,2.5,2,1), mgp=c(1.6,.6,0))
plot(ACF,  type="h", xlab="LAG", ylim=c(-.4,.8), axes=FALSE)  
abline(h=0)
axis(1, seq(0,50,by=12))
axis(2)
box()
plot(PACF, type="h", xlab="LAG", ylim=c(-.4,.8), axes=FALSE)  
abline(h=0)
axis(1, seq(0,50,by=12))
axis(2)
box()

dev.print(pdf,paste(c(path_fig,"sarmaacf.pdf"),collapse="/"))
quartz_off_screen 
                2 

SARIMA

http://www.personal.psu.edu/asb17/old/sta4853/files/sta4853-17print.pdf

autoplot(prodn)

dev.print(pdf,paste(c(path_fig,"prodn.pdf"),collapse="/"))
quartz_off_screen 
                2 
autoplot(diff(prodn))

dev.print(pdf,paste(c(path_fig,"diffprodn.pdf"),collapse="/"))
quartz_off_screen 
                2 
monthplot(diff(prodn))

dev.print(pdf,paste(c(path_fig,"months_diffprodn.pdf"),collapse="/"))
quartz_off_screen 
                2 
ggAcf(diff(prodn),lag.max = 40)

dev.print(pdf,paste(c(path_fig,"acf_diffprodn.pdf"),collapse="/"))
quartz_off_screen 
                2 
ggPacf(diff(prodn),lag.max = 40)

dev.print(pdf,paste(c(path_fig,"pacf_diffprodn.pdf"),collapse="/"))
quartz_off_screen 
                2 
LS0tCnRpdGxlOiAiQ2hhcHRlciA2OiBCdWlsZGluZyBBUklNQSBhbmQgU0FSSU1BIG1vZGVscyIKb3V0cHV0OiBodG1sX25vdGVib29rCi0tLQoKYGBge3J9CmxpYnJhcnkoZ2dmb3J0aWZ5KQpsaWJyYXJ5KGZvcmVjYXN0KQpwYXRoX2ZpZyA9ICJ+L0Ryb3Bib3gvRXZyeS9UaW1lU2VyaWVzL0NvdXJzL0ZpZ3VyZXMiCmBgYAoKIyBBbiBleGFtcGxlIG9mIEFSSU1BIG1vZGVsIApXZSBjb25zaWRlciB0aGUgYW5hbHlzaXMgb2YgcXVhcnRlcmx5IFUuUy4gR05QIGZyb20KMTk0NygxKSB0byAyMDAyKDMpLCBuID0gMjIzIG9ic2VydmF0aW9ucy4gVGhlIGRhdGEgYXJlIHJlYWwgVS5TLiBncm9zcwpuYXRpb25hbCBwcm9kdWN0IGluIGJpbGxpb25zIG9mIGNoYWluZWQgMTk5NiBkb2xsYXJzIGFuZCBoYXZlIGJlZW4gc2Vhc29uYWxseSBhZGp1c3RlZC4gVGhlIGRhdGEgd2VyZSBvYnRhaW5lZCBmcm9tIHRoZSBGZWRlcmFsIFJlc2VydmUgQmFuayBvZiBTdC4gTG91aXMgKGh0dHA6Ly9yZXNlYXJjaC5zdGxvdWlzZmVkLm9yZy8pLiBUaGlzIGV4YW1wbGUgaXMgZnJvbSBQZXRlciBKIEJyb2Nrd2VsbCBhbmQgUmljaGFyZCBBIERhdmlzLCBUaW1lIHNlcmllcyA6IHRoZW9yeSBhbmQgbWV0aG9kcywgU3ByaW5nZXIgU2NpZW5jZSAmIEJ1c2luZXNzIE1lZGlhLCAyMDEzLgoKSGVyZSBhcmUgdGhlIHN0ZXBzIHRvIGJ1aWxkIGFuIEFSSU1BIG1vZGVsIGZvciB0aGUgYGducGAgc2VyaWVzIG9mIGBhc3RzYWAgbGlicmFyeS4KYGBge3J9CmxpYnJhcnkoYXN0c2EpCmBgYAoKIyMgMS4gVGltZSBwbG90IG9mIGBnbnBgCmBgYHtyfQphdXRvcGxvdChnbnApCmBgYApUaGVyZSBpcyBhIGNsZWFyIHRyZW5kLCB0aGUgdmFyaWFiaWxpdHkgaXMgcGVyaGFwcyBpbmNyZWFzaW5nIHdpdGggdGhlIG1lYW4uIFdlJ2xsIHRoZSBkaWZmZXJlbnRpYXRpb24gb24gYGxvZyhnbnApYCBhbmQgYGducGAKYGBge3J9CmdnQWNmKGducCkKYGBgClRoZSBBQ0ZzIGlzLCBhcyBleHBlY3RlZCwgKHRvbykgc2xvd2x5IGRlY3JlYXNpbmcuCgojIyAyLiBUcmFuc2Zvcm1hdGlvbiBvZiB0aGUgZGF0YQpgYGB7cn0KYXV0b3Bsb3QobG9nKGducCkpCmdnQWNmKGxvZyhnbnApKQpgYGAKCiMjIDMuIERpZmZlcmVudGlhdGlvbgoKIyMjIDMuMSBGb3IgYGducGAKYGBge3J9CmF1dG9wbG90KGRpZmYoZ25wKSkKZ2dBY2YoZGlmZihnbnApKQpgYGAKClRoaXMgaXMgc3RpbGwgbm90IHN0YXRpb25hcnk6IG5laXRoZXIgYXQgb3JkZXIgMSAobWVhbiksIG5vciBhdCBvcmRlciAyICh2YXJpYW5jZSkuIExldCdzIHRyeSBoaWdoZXIgb3JkZXJzLgpgYGB7cn0KYXV0b3Bsb3QoZGlmZihnbnAsMikpCmdnQWNmKGRpZmYoZ25wLDIpKQpgYGAKCiMjIyAzLjIgRm9yIGBsb2coZ25wKWAKYGBge3J9CmF1dG9wbG90KGRpZmYobG9nKGducCkpKQpnZ0FjZihkaWZmKGxvZyhnbnApKSkKYGBgClRoaXMgbG9va3MgcmVhc29uYWJpbGl0eSBzdGF0aW9ubmFyeSBhdCBvcmRlcnMgMSBhbmQgMi4gV2Ugbm93IGNvbnN0cnVjdCBhbiBBUk1BIG1vZGVsIGZvciBgZGlmZihsb2coZ25wKSlgCgojIyA0LiBBQ0ZzIGFuZCBQQUNGcyBmb3IgYGRpZmYobG9nKGducCkpYAoKYGBge3J9CmdnQWNmKGRpZmYobG9nKGducCkpKQpnZ1BhY2YoZGlmZihsb2coZ25wKSkpCmBgYApUaGUgQUNGIGNhbGxzIGZvciBhbiBNQSgyKSwgd2hpbGUgUEFDRiBjYWxscyBmb3IgYW4gQVIoMSkuCldlJ2xsIHRoZW4gZXN0aW1hdGUgaW4gQVJJTUEoMSwxLDApIGFuZCBBUklNQSgwLDEsMikgbW9kZWxzIGZvciBgbG9nKGducClgLgoKIyMgNS4gTW9kZWwgc2VsZWN0aW9uCgoqKkNhdXRpb24qKjogdGhlIGRlZmF1bHQgZm9yIHRoZSBvcHRpb24gYGluY2x1ZGUuZHJpZnRgIGluIHRoZSBgQXJpbWFgIGZ1bmN0aW9uIG9mIHRoZSBgZm9yZWNhc3RgIHBhY2thZ2UgaXMgRkFMU0UsIGJ1dCBpbiBvdXIgY2FzZSB0aGUgYGRpZmYobG9nKGducCkpYCBpcyBjbGVhcmx5IG5vdCBjZW50ZXJlZApgYGB7cn0KYXV0b3Bsb3QoZGlmZihsb2coZ25wKSkpCmBgYAoKYGBge3J9CmFyaW1hXzExMCA9IEFyaW1hKGxvZyhnbnApLG9yZGVyID0gYygxLDEsMCksaW5jbHVkZS5jb25zdGFudD1UUlVFKQphcmltYV8wMTIgPSBBcmltYShsb2coZ25wKSxvcmRlciA9IGMoMCwxLDIpLGluY2x1ZGUuY29uc3RhbnQ9VFJVRSkKY3JpdGVyaWEgPSBtYXRyaXgoYyhhcmltYV8xMTAkYWljYyxhcmltYV8wMTIkYWljYyxhcmltYV8xMTAkYmljLGFyaW1hXzAxMiRiaWMpLG5yb3c9MixieXJvdyA9IFRSVUUpCmNvbG5hbWVzKGNyaXRlcmlhKSA9IGMoImFyaW1hXzExMCIsImFyaW1hXzAxMiIpCnJvd25hbWVzKGNyaXRlcmlhKSA9IGMoImFpY2MiLCJiaWMiKQpjcml0ZXJpYQpgYGAKCiMjIDYuIERpYWdub3NpcyBvbiByZXNpZHVhbHMKYGBge3J9CmNoZWNrID0gY2hlY2tyZXNpZHVhbHMoYXJpbWFfMTEwICkKYGBgCmBgYHtyfQphdXRvLmFyaW1hKGxvZyhnbnApKQpgYGAKCiMgQW4gZXhhbXBsZSBvZiBTQVJJTUEgbW9kZWwKYGBge3J9CgpgYGAKCiMjIEEgcHVyZSBwdXJlIHNlYXNvbmFsICRBUk1BKDAsMSlfUyQgbW9kZWwKYGBge3J9CnNldC5zZWVkKDY2NikKb21lZ2EgPXJub3JtKG49MTAwMCkKWCA9IHRzKG9tZWdhWzEzOjEwMDBdIC0gMC41ICogb21lZ2FbMTo5ODhdLCBzdGFydCA9IGMoMTkwMCwxKSwgZnJlcXVlbmN5ID0gMTIpCmF1dG9wbG90KFgpCmdnQWNmKFgpCmRldi5wcmludChwZGYscGFzdGUoYyhwYXRoX2ZpZywiYWNmX3B1cmVfc2Vhc29uX01BMS5wZGYiKSxjb2xsYXBzZT0iLyIpKQpnZ1BhY2YoWCkKZGV2LnByaW50KHBkZixwYXN0ZShjKHBhdGhfZmlnLCJwYWNmX3B1cmVfc2Vhc29uX01BMS5wZGYiKSxjb2xsYXBzZT0iLyIpKQpgYGAKYGBge3J9CnBoaSA9IGMocmVwKDAsMTEpLC44KQpBQ0YgPSBBUk1BYWNmKGFyPXBoaSwgbWE9LS41LCA1MClbLTFdICAgICAjIFstMV0gcmVtb3ZlcyAwIGxhZwpQQUNGID0gQVJNQWFjZihhcj1waGksIG1hPS0uNSwgNTAsIHBhY2Y9VFJVRSkKcGFyKG1mcm93PWMoMSwyKSwgbWFyPWMoMi41LDIuNSwyLDEpLCBtZ3A9YygxLjYsLjYsMCkpCnBsb3QoQUNGLCAgdHlwZT0iaCIsIHhsYWI9IkxBRyIsIHlsaW09YygtLjQsLjgpLCBheGVzPUZBTFNFKSAgCmFibGluZShoPTApCmF4aXMoMSwgc2VxKDAsNTAsYnk9MTIpKQpheGlzKDIpCmJveCgpCnBsb3QoUEFDRiwgdHlwZT0iaCIsIHhsYWI9IkxBRyIsIHlsaW09YygtLjQsLjgpLCBheGVzPUZBTFNFKSAgCmFibGluZShoPTApCmF4aXMoMSwgc2VxKDAsNTAsYnk9MTIpKQpheGlzKDIpCmJveCgpCmRldi5wcmludChwZGYscGFzdGUoYyhwYXRoX2ZpZywic2FybWFhY2YucGRmIiksY29sbGFwc2U9Ii8iKSkKCgpgYGAKCiMjIFNBUklNQQpodHRwOi8vd3d3LnBlcnNvbmFsLnBzdS5lZHUvYXNiMTcvb2xkL3N0YTQ4NTMvZmlsZXMvc3RhNDg1My0xN3ByaW50LnBkZgoKYGBge3J9CmF1dG9wbG90KHByb2RuKQpkZXYucHJpbnQocGRmLHBhc3RlKGMocGF0aF9maWcsInByb2RuLnBkZiIpLGNvbGxhcHNlPSIvIikpCmBgYApgYGB7cn0KYXV0b3Bsb3QoZGlmZihwcm9kbikpCmRldi5wcmludChwZGYscGFzdGUoYyhwYXRoX2ZpZywiZGlmZnByb2RuLnBkZiIpLGNvbGxhcHNlPSIvIikpCmBgYAoKCmBgYHtyfQptb250aHBsb3QoZGlmZihwcm9kbikpCmRldi5wcmludChwZGYscGFzdGUoYyhwYXRoX2ZpZywibW9udGhzX2RpZmZwcm9kbi5wZGYiKSxjb2xsYXBzZT0iLyIpKQpgYGAKCgoKCgpgYGB7cn0KZ2dBY2YoZGlmZihwcm9kbiksbGFnLm1heCA9IDQwKQpkZXYucHJpbnQocGRmLHBhc3RlKGMocGF0aF9maWcsImFjZl9kaWZmcHJvZG4ucGRmIiksY29sbGFwc2U9Ii8iKSkKCmBgYAoKCmBgYHtyfQpnZ1BhY2YoZGlmZihwcm9kbiksbGFnLm1heCA9IDQwKQpkZXYucHJpbnQocGRmLHBhc3RlKGMocGF0aF9maWcsInBhY2ZfZGlmZnByb2RuLnBkZiIpLGNvbGxhcHNlPSIvIikpCmBgYAoKCg==