ARIMA
thomas petersen
9/7/2017
Gennemgangen af tidsrækkeanalyse bygger meget på praktisk anvendelighed (dvs. vi vil gerne kunne forudsige kursudviklingen), vi vil springe let hen over teorien der kan være tung og er meget omfattende. Nedenstående links giver dog en indføring i den teoretiske del, som vi her ikke berører.
http://ucanalytics.com/blogs/arima-models-manufacturing-case-study-example-part-3/
Her er en gennemgang af forskellige typer af tidsrækker man kan opleve.
https://people.duke.edu/~rnau/411arim.htm#arima010
Video om ARIMA https://youtu.be/Aw77aMLj9uM
ARIMA(0,0,0)
Nedenfor har vi aktiekurser for 50 dage for en fiktiv aktie, vi vil nu undersøge om disse kan bruges til at forudsige noget om fremtidige aktiekurser. For at gøre dette, skal man enten importere Excelfilen, eller copy paste data fra rammen nedenfor.
Data nedenfor kan kopieres direkte ind i R. Navngiv data ARIMA1 ved at skrive følgende tekst foran det du kopierer.
ARIMA1 <-
structure(c(21.64, 19.62, 14.7, 20.52, 19.94, 20.1, 19.46, 20.76,
20.64, 15.66, 18.76, 22.24, 22.94, 16.56, 18.44, 17.98, 20.6,
19.86, 18.28, 18.1, 18.54, 20.86, 19.4, 23.86, 18.44, 20.48,
21.92, 19.42, 18.96, 22.52, 21.28, 19.64, 21.38, 18.26, 16.48,
17.92, 20.72, 20.14, 18.58, 21.34, 21.12, 21.96, 19.74, 20.02,
16.3, 23.12, 21.26, 18.52, 24.58, 21.26), .Dim = c(50L, 1L), .Dimnames = list(
NULL, "Dagskurs"), .Tsp = c(1, 50, 1), class = "ts")Vi kan nu plotte vore data i R.
plot.ts(ARIMA1, xlab='Tid', ylab = 'Kursdata')
Det er svært at se nogen tydelig udvikling i kursen.
Vi benytter auto.arima til at undersøge om der er en systematik i tidsserien, for at bruge denne funktion skal vi hente og loade pakken forecast med hhv. install.packages(“forecast”) kun første gang og require(“forecast”) eller library(“forecast”)
## Loading required package: forecastauto.arima(ARIMA1)## Series: ARIMA1
## ARIMA(0,0,0) with non-zero mean
##
## Coefficients:
## mean
## 19.8964
## s.e. 0.2872
##
## sigma^2 estimated as 4.209: log likelihood=-106.37
## AIC=216.74 AICc=217 BIC=220.57Funktionen auto.arima i R finder automatisk den ARIMA model der passer bedst på observationerne.
Output ARIMA(0,0,0) with non-zero mean, fortæller os at data er ligesom hvid støj. Den bedste forudsigelse af aktieprisen vi kan komme med er gennemsnittet af alle kurserne. Vi kan altså ikke forudsige prisen vha. vore fine værktøjer.
Akaike Information Criterion (AIC) , og Bayesian Information Criterion (BIC) benyttes til at vælge ARIMA modellen med mindst AIC og BIC værdier. auto.arima finder den bedste model automatisk.
Her er ligningen for aktiekursen, den bedste forudsigelse af den fremtidige kurs er den gennemsnitlige kurs der tidligere er observeret.
\[\hat{Y_t}=19.90\]
ARIMA(1,0,0) eller AR(1) autoregression
Vi ser nu på 50 dagskurser for en anden aktie, kopier data ind i R og navngiv med ARIMA2:
structure(c(97.957, 73.334, 83.723, 92.322, 106.694, 94.278,
118.855, 107.169, 116.21, 124.179, 105.945, 100.359, 131.513,
103.029, 104.339, 103.513, 108.638, 97.032, 100.426, 91.356,
110.501, 116.041, 106.145, 96.257, 101.983, 124.107, 112.892,
121.101, 103.067, 98.862, 94.009, 115.093, 103.553, 92.895, 103.087,
89.622, 91.413, 116.587, 119.88, 111.579, 100.377, 108.997, 104.367,
82.48, 79.99, 85.062, 102.597, 98.046, 101.043, 109.051), .Dim = c(50L,
1L), .Dimnames = list(NULL, "Dagskurs"), .Tsp = c(1, 50, 1), class = "ts")ts.plot(ARIMA2)
auto.arima(ARIMA2)## Series: ARIMA2
## ARIMA(1,0,0) with non-zero mean
##
## Coefficients:
## ar1 mean
## 0.3664 103.2372
## s.e. 0.1301 2.4492
##
## sigma^2 estimated as 128.3: log likelihood=-191.36
## AIC=388.73 AICc=389.25 BIC=394.46Her afslører auto.arima 1. ordens autoregression dvs.
Modellen kan skrives som.
\[\hat{Y_t}=c + \phi Y_{t-1}\] Vi forestiller os eksemplet hvor den sande middelværdi for tidsrækken er \(\mu=100\) og \(\phi=0.5\) for en ARIMA(1,0,0) model. Så kan vi beregne \(c=(1-\phi)\cdot \mu=(1-0.5)\cdot 100=50\) er middelværdien estimeret ved den gennemsnitlige kurs. Ligningen for modellen kan så skrives som:
\[\hat{Y_t}=c + \phi Y_{t-1} \Leftrightarrow \hat{Y_t}=50 + 0.5 Y_{t-1}\]
\(\phi\) fortæller, hvis kursen dagen før var 80 gennemsnitskursen er 100, vil kursen imorgen ifølge modellen være forudsagt som \(50+0.5\cdot 80=90\) dagen efter vil den så være \(50+0.5\cdot 90=95\). De forudsagte værdier konvergerer (nærmer sig) mod \(\mu\). AR i ARIMA, står for autoregression, selv-regression mod middelværdien, i eksemplet så vi hvordan værdien nærmer sig 100, hvis vi forudsiger flere dages kurser kan vi se dette.
\(\phi\) må kun antage værdier mellem og ikke lig med -1 og 1, hvilket betyder den er stationær, altså nærmer sig den sande middelværdi \(\mu\).
Hvad vil der ske hvis \(\mu=100\) og \(\phi=-0.5\) for en ARIMA(1,0,0) model (husk \(c=(1-\phi)\cdot\mu\) når man skal bestemme modellen)?
Vi kan nu i R forudsige aktiekursen 12 perioder frem med predict:
predict(auto.arima(ARIMA2), n.ahead = 12)$pred## Time Series:
## Start = 51
## End = 62
## Frequency = 1
## [1] 105.3673 104.0177 103.5232 103.3420 103.2756 103.2513 103.2424
## [8] 103.2391 103.2379 103.2375 103.2373 103.2373ARIMA(0,1,0) eller I(1) Random Walk with a drift
Hvis en serie er ikke-stationær, er den simpleste model en random walk:
\[\hat{Y_t}-Y_{t-1}=\mu\Leftrightarrow \hat{Y_t}=\mu-Y_{t-1}\] Dette betyder at Y stiger konstant med \(\mu\) i hver periode. Drift betyder at tidsrækken stiger konstant.
Kald den næste vektor for ARIMA3:
structure(c(100, 110.695, 116.16, 123.983, 124.182, 118.432,
114.157, 116.126, 126.584, 139.104, 149.709, 166.794, 172.85,
168.264, 184.061, 188.304, 196.156, 203.239, 217.565, 226.016,
236.986, 250.127, 254.33, 255.305, 280.172, 295.641, 323.152,
337.009, 344.1, 371.613, 383.138, 390.722, 390.121, 392.329,
408.905, 423.546, 411.232, 436.033, 441.104, 433.005, 436.042,
455.257, 455.92, 473.371, 475.472, 490.525, 492.407, 500.438,
506.987, 504.54, 495.746), .Dim = c(51L, 1L), .Dimnames = list(
NULL, "Dagskurs"), .Tsp = c(1, 51, 1), class = "ts")ts.plot(ARIMA3)
auto.arima(ARIMA3)## Series: ARIMA3
## ARIMA(0,1,0) with drift
##
## Coefficients:
## drift
## 7.9149
## s.e. 1.2790
##
## sigma^2 estimated as 83.46: log likelihood=-181.05
## AIC=366.1 AICc=366.36 BIC=369.93Modellen ovenfor kan skrives som: \[\hat{Y_t}-Y_{t-1}=\mu\Leftrightarrow \hat{Y_t}-Y_{t-1}=7.9\] Vi indsætter drift i stedet for \(\mu\), tolningen er at modellen forudsiger at aktiekursen stiger med 7.9 fra periode til periode.
Forstiller man sig en ARIMA(0,1,0) med drift 10 og en kurs på tidspunkt t-1 på 120, vil vi forudsige en kurs på 130 ved tid t og 140 ved tid t+1 osv.
Hvis vi har en ren random walk model uden drift dvs. med \(\mu=0\) ARIMA(0,1,0) for en aktiekurs , forventer vi at kursen til tid t vil være den samme som til tid t-1. Denne kan skrives som:
\[\hat{Y_t}-Y_{t-1}=0\]
ARIMA(0,0,1) eller MA(1) Moving average
Vi kan i stedet for at bruge tidligere aktiekurser til at forudsige aktiekursen i stedet benytte tidligere målefejl residualer til at forudsige kursen.
Modellen kan skrives som:
\[\hat{Y_t}=\mu+\theta_1 e_{t-1}\] Hvis vi forestiller os \(\mu=50\) \(\theta_1=0.5\) kursen til tid t-1 var 120 forudsigelsen til tid t-1 var 100, så målefejlen residualen til tid t-1 er \(e_{t-1}\) er faktisk kurs minus forudsagt kurs altså 120-100=20. Nu kan vi forudsige kursen til tid t som: \[\hat{Y_t}=\mu+\theta_1 e_{t-1}\Leftrightarrow \hat{Y_t}=50+0.5\cdot20=60\]
Data der kan indlæses i R, benævn data som ARIMA4:
structure(c(104.035, 84.696, 66.622, 75.078, 90.512, 80.213,
95.646, 120.728, 107.473, 98.007, 105.692, 91.47, 81.564, 90.325,
112.94, 128.034, 118.113, 108.337, 92.803, 85.981, 91.068, 90.237,
94.939, 87.966, 93.857, 110.179, 99.252, 84.86, 80.571, 74.44,
107.43, 112.259, 97.252, 96.04, 93.651, 113.215, 116.992, 100.672,
93.411, 104.566, 115.626, 93.047, 91.644, 120.28, 109.674, 95.901,
102.858, 105.829, 99.581, 112.855), .Dim = c(50L, 1L), .Dimnames = list(
NULL, "Dagskurs"), .Tsp = c(1, 50, 1), class = "ts")ts.plot(ARIMA4)
auto.arima(ARIMA4)## Series: ARIMA4
## ARIMA(0,0,1) with non-zero mean
##
## Coefficients:
## ma1 mean
## 0.9053 99.0176
## s.e. 0.0664 2.5588
##
## sigma^2 estimated as 95.68: log likelihood=-184.81
## AIC=375.62 AICc=376.14 BIC=381.35Vi kan nu forudsige aktiekursen 12 perioder frem med predict:
predict(auto.arima(ARIMA4), n.ahead = 12)$pred## Time Series:
## Start = 51
## End = 62
## Frequency = 1
## [1] 113.64960 99.01756 99.01756 99.01756 99.01756 99.01756 99.01756
## [8] 99.01756 99.01756 99.01756 99.01756 99.01756Hvorfor svarer den forudsagte værdi til mean i en ren ARIMA(0,0,1) eller MA(1) model? (Vink hvad er definitionen på en residual)
Plots med forskellige modeller
## Series: ap1
## ARIMA(0,1,1)
##
## Coefficients:
## ma1
## -0.3783
## s.e. 0.1455
##
## sigma^2 estimated as 111.4: log likelihood=-188.35
## AIC=380.7 AICc=380.96 BIC=384.53Kursen svinger omkring middelværdien.
## Series: ap2
## ARIMA(0,1,1) with drift
##
## Coefficients:
## ma1 drift
## 0.6162 -5.406
## s.e. 0.1303 1.768
##
## sigma^2 estimated as 63.27: log likelihood=-173.85
## AIC=353.7 AICc=354.22 BIC=359.43
## Series: ap3
## ARIMA(2,0,1) with non-zero mean
##
## Coefficients:
## ar1 ar2 ma1 mean
## -0.3225 0.1403 0.8435 98.1645
## s.e. 0.1128 0.0924 0.0837 1.0982
##
## sigma^2 estimated as 101.6: log likelihood=-744.14
## AIC=1498.28 AICc=1498.59 BIC=1514.77Vi kan også grafisk vise hvordan kursen vil udvikle sig med 80% og 95% konfidensbælter.
forecast(auto.arima(ap3))## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 201 98.70266 85.78395 111.6214 78.94520 118.4601
## 202 97.23784 82.67080 111.8049 74.95947 119.5162
## 203 98.53885 83.96741 113.1103 76.25376 120.8239
## 204 97.91370 83.30376 112.5236 75.56972 120.2577
## 205 98.29789 83.68269 112.9131 75.94586 120.6499
## 206 98.08626 83.46846 112.7041 75.73027 120.4422
## 207 98.20842 83.58992 112.8269 75.85135 120.5655
## 208 98.13933 83.52058 112.7581 75.78188 120.4968
## 209 98.17875 83.55993 112.7976 75.82119 120.5363
## 210 98.15634 83.53749 112.7752 75.79873 120.5139plot(forecast(auto.arima(ap3)))
ARIMA af højere orden
Arima modeller kan afhænge af flere tidligere perioder, fx kan ligningen for ARIMA(2,0,0) eller AR(2), opskrives som:
\[\hat{Y_t}=c + \phi Y_{t-1}+ \phi_2 Y_{t-2}\] Modellen afhænger altså af 2 tidligere perioder (lags) og ikke en. Man betegner dette som en model med lag 2.
Arima modeller kan indeholde flere forskellige elementer med lag som fx. ARIMA(0,2,1).
ARIMA og sæsonalitet
Hvis fx. en aktie handles lavere om fredagen kan ARIMA modellerne korrigere for dette ved sæsonkorrektion. I sæsonkorrigerede modeller vises dette som en ekstra vektor med 3 tal for hhv. sæsonkorrigeret AR eller SAR, sæsonkorrigeret I eller SI og sæsonkorrigeret MA eller SMA. En model som ARIMA(1,0,0)(1,0,0) har altså udover AR også en sæsonkomponent.
ARIMA eksempler
Traktorer
Hent følgende data for traktor salg, med følgende kommandoer i R.
data = read.csv('http://ucanalytics.com/blogs/wp-content/uploads/2015/06/Tractor-Sales.csv')
data = ts(data[,2],start = c(2003,1),frequency = 12)Vi ser salget er voksende over tid, der er ligeledes en sæsonkomponent.
plot(data, xlab='Years', ylab = 'Tractor Sales')
Differens tranformer data for at generere stationære data mht. middel (fjern trend)
plot(diff(data),ylab='Differenced Tractor Sales')
log transformer data for at sikre stationaritet mht. varians.
plot(log10(data),ylab='Log (Tractor Sales)')
Eventuel Differens og log transformation af data for at sikre stationaritet både mht. middel og varians.
plot(diff(log10(data)),ylab='Differenced Log (Tractor Sales)')
Find bedste model med auto.arima, når der er stationaritet.
Akaike Information Criterion (AIC) , og Bayesian Information Criterion (BIC), vælg ARIMA modellen med mindst AIC and BIC værdier. auto.arima finder den bedste model automatisk.
require(forecast)
ARIMAfit = auto.arima(log10(data), approximation=FALSE,trace=FALSE)
ARIMAfit## Series: log10(data)
## ARIMA(0,1,1)(0,1,1)[12]
##
## Coefficients:
## ma1 sma1
## -0.4047 -0.5529
## s.e. 0.0885 0.0734
##
## sigma^2 estimated as 0.0002571: log likelihood=354.4
## AIC=-702.79 AICc=-702.6 BIC=-694.17Nu kan vi forudsige kommende traktor salg med modellen
par(mfrow = c(1,1))
pred = predict(ARIMAfit, n.ahead = 36)
salg <- 10^pred$pred
salg## Jan Feb Mar Apr May Jun Jul
## 2015 567.7645 566.4765 670.8226 758.9138 855.9482 817.2827 938.7239
## 2016 625.2464 623.8280 738.7384 835.7481 942.6065 900.0265 1033.7626
## 2017 688.5479 686.9859 813.5300 920.3613 1038.0383 991.1474 1138.4233
## Aug Sep Oct Nov Dec
## 2015 934.5120 703.5005 626.9879 571.9988 668.5363
## 2016 1029.1243 774.7246 690.4657 629.9094 736.2206
## 2017 1133.3154 853.1596 760.3701 693.6830 810.7573plot(data,type='l',xlim=c(2003,2018),ylim=c(1,1600),xlab = 'Year',ylab = 'Tractor Salg')
lines(10^(pred$pred),col='blue')
lines(10^(pred$pred+2*pred$se),col='orange')
lines(10^(pred$pred-2*pred$se),col='orange')
Detail debet card forbrug på Island (millioner ISK).
#Hent fpp pakken og load den
plot(debitcards)
dldebitcards <- diff(log10(debitcards))
plot(dldebitcards,ylab="Differenced Log (debitcards)")
require(forecast)
ARIMAfit = auto.arima(log10(debitcards), approximation=FALSE,trace=FALSE)
ARIMAfit## Series: log10(debitcards)
## ARIMA(0,1,1)(1,0,0)[12]
##
## Coefficients:
## ma1 sar1
## -0.6843 0.8747
## s.e. 0.0574 0.0328
##
## sigma^2 estimated as 0.0006647: log likelihood=339.09
## AIC=-672.18 AICc=-672.02 BIC=-663.05Nu kan vi forudsige kommende debetkort omsætning med modellen
par(mfrow = c(1,1))
pred = predict(ARIMAfit, n.ahead = 36)
plot(debitcards,type='l',xlim=c(2000,2016),ylim=c(1,40000),xlab = 'Year',ylab = 'Debetcard usage')
lines(10^(pred$pred),col='blue')
lines(10^(pred$pred+2*pred$se),col='orange')
lines(10^(pred$pred-2*pred$se),col='orange')
Forudsagt brug af debetkort bliver:
10^(pred$pred)## Jan Feb Mar Apr May Jun Jul
## 2013 19742.32 19660.50 20801.63 21169.26 23415.54 23427.33 24136.87
## 2014 20818.53 20743.04 21792.45 22128.98 24169.86 24180.50 24819.93
## 2015 21807.94 21738.75 22697.80 23004.11 24849.60 24859.18 25433.27
## Aug Sep Oct Nov Dec
## 2013 25522.95 21355.03 22921.36 21963.15 27088.52
## 2014 26062.32 22298.76 23723.06 22853.24 27455.49
## 2015 26543.48 23158.42 24447.30 23661.38 27780.58Forecast Aktiekurser
Man kan hente online aktiekurser med quantmod pakken installer denne med install.packages og library kommandoerne. Vi henter nedenfor Google lukkekurs til dato.
getSymbols("GOOG",from = "2017-01-01", to = Sys.Date(),getSymbols.warning4.0=FALSE)## [1] "GOOG"plot(GOOG[,6],main = "Google adj. close")
agoog <- auto.arima(GOOG[,6])
agoog## Series: GOOG[, 6]
## ARIMA(0,1,0) with drift
##
## Coefficients:
## drift
## 1.0039
## s.e. 0.6260
##
## sigma^2 estimated as 75.63: log likelihood=-687.21
## AIC=1378.43 AICc=1378.49 BIC=1384.94fagoog <- forecast(agoog)
fagoog## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 194 979.8939 968.7488 991.0391 962.8489 996.939
## 195 980.8978 965.1362 996.6594 956.7925 1005.003
## 196 981.9017 962.5978 1001.2057 952.3789 1011.425
## 197 982.9056 960.6153 1005.1959 948.8156 1016.996
## 198 983.9095 958.9882 1008.8309 945.7957 1022.023
## 199 984.9135 957.6135 1012.2134 943.1618 1026.665
## 200 985.9174 956.4301 1015.4046 940.8204 1031.014
## 201 986.9213 955.3980 1018.4445 938.7106 1035.132
## 202 987.9252 954.4897 1021.3606 936.7901 1039.060
## 203 988.9291 953.6850 1024.1731 935.0280 1042.830plot(fagoog,main = "Google adj. close")
require("quantmod")
getSymbols("GS",from = "2017-01-01", to = Sys.Date(),getSymbols.warning4.0=FALSE)## [1] "GS"plot(GS[,6],main = "Goldman Sachs adj. close")
ags <- auto.arima(GS[,6])
ags## Series: GS[, 6]
## ARIMA(0,1,0)
##
## sigma^2 estimated as 8.387: log likelihood=-476.6
## AIC=955.19 AICc=955.22 BIC=958.45fags <- forecast(ags)
fags## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 194 246.02 242.3085 249.7315 240.3438 251.6962
## 195 246.02 240.7712 251.2688 237.9926 254.0474
## 196 246.02 239.5915 252.4485 236.1885 255.8515
## 197 246.02 238.5971 253.4429 234.6676 257.3724
## 198 246.02 237.7209 254.3191 233.3276 258.7124
## 199 246.02 236.9288 255.1112 232.1162 259.9238
## 200 246.02 236.2004 255.8396 231.0022 261.0378
## 201 246.02 235.5224 256.5176 229.9653 262.0747
## 202 246.02 234.8856 257.1544 228.9914 263.0486
## 203 246.02 234.2833 257.7567 228.0703 263.9697plot(fags,main = "Goldman Sachs adj. close")
require("quantmod")
getSymbols("DANSKE.CO",from = "2017-01-01", to = Sys.Date(),getSymbols.warning4.0=FALSE)## [1] "DANSKE.CO"plot(DANSKE.CO[,6],main = "Danske Bank adj. close")
addb <- auto.arima(DANSKE.CO[,6])
addb## Series: DANSKE.CO[, 6]
## ARIMA(0,1,0)
##
## sigma^2 estimated as 4.967: log likelihood=-426.31
## AIC=854.62 AICc=854.64 BIC=857.88faddb <- forecast(addb)
faddb## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 194 250.8 247.9437 253.6563 246.4317 255.1683
## 195 250.8 246.7606 254.8394 244.6223 256.9777
## 196 250.8 245.8528 255.7472 243.2339 258.3661
## 197 250.8 245.0875 256.5125 242.0634 259.5366
## 198 250.8 244.4132 257.1868 241.0322 260.5678
## 199 250.8 243.8036 257.7964 240.0999 261.5001
## 200 250.8 243.2430 258.3570 239.2426 262.3574
## 201 250.8 242.7213 258.8787 238.4446 263.1554
## 202 250.8 242.2312 259.3688 237.6952 263.9048
## 203 250.8 241.7677 259.8323 236.9863 264.6137plot(faddb,main = "Danske Bank adj. close")