Dynaaminen regressiomalli

Kun haluamme sovittaa regressiomallin aikasarjadataan, niin on tärkeää huomioida aikasarjoille tyypillinen epästationaarinen rakenne, eli rakenne, jossa aikasarjan aiempi arvo ennustaa sen tulevaa arvoa. Alla on esimerkki epästationaarisesta aikasarjasta, jossa jälkimmäisen kuukauden arvo on lähes aina aiemman kuukauden arvoa suurempi.

Lähde: http://www.uow.edu.au/student/qualities/statlit/module3/9time/index.html

Jos tällaiseen dataan sijoittaa regressiomallin suoraan, niin sekä regressiokertoimet että keskivirheet ovat vääristyneitä. Tällöin on kyse ns. Spurious Regression:sta, jolloin saisimme pienet p-arvot laittamalla muuttujaksi esim. maapallon keskilämpötilan kehityksen samana ajankohtana. Spurious regression siis altistaa mallin vääränlaisille tulkinnoille. Aikasarjadataa käsitellessä on siis tärkeää tutkia ensin, onko data epästationaarisessa muodossa vai ei. Datan saa stationaariseen muotoon differensoimalla sen tarpeeksi monta kertaa. Yleensä 1-2 kertaa riittää. Jos datan muodostaneen AR(1)-termin kerroin on yli yksi, niin data voi tarvita kolme differensointia. Alla on esimerkki datasta, joka vaatii kolme differensointia:
R-koodi:

z <-rnorm(100,1,0)
g <- rnorm(100,mean=0,sd=1)
for (t in 2:100) z[t]=1.1*z[t-1]+g[t]

Dynaaminen regressiomalli

Käsittelen nyt tilannetta, jossa regressiomallin virhetermi seuraa ARIMA-mallia, eli regression tehokkuusehdot eivät täyty, koska sen virhetermit korreloivat keskenään. Tällöin emme hyödynnä kaikkea informaatiota mallissamme ja lisäksi keskivirheet ovat vääristyneitä. On siis tärkeää hyödyntää ARIMA-mallia yhdessä regressiomallin kanssa. Muodostan ensin datan, jossa y on riippuvainen x:sta sekä ARIMA-mallia seuraavasta virhetermistä.
R-koodi

y  <- rnorm(100,mean=0,sd=1)
g <- rnorm(100,mean=0,sd=1)
x <- rnorm(100,mean=0,sd=5)
e <- rnorm(100,mean=0,sd=1)
for (t in 2:100) e[t]=e[t-1]+g[t] # Muodostan yksikköjuurellisen (seuraa ARIMA-mallia) virhetermin
for (t in 2:100) y[t]=x[t]+e[t]
library(tseries)
adf.test(y) # Dickey-Fuller -testin mukaan y on stationaarinen.

Seuraavaksi muodostan regressiomallin, jolla saamme x:n kertoimeksi 1.040 (Keskivirhe on 0.038). Testaan sitten virhetermien autokorrelaation Breusch-Godfrey -testillä, jotta näen, ovatko mallin tulokset luotettavia. Testin p-arvo on pienempi kuin 0.05, mikä viittaa virhetermeissä esiintyvään autokorrelaatioon. Alla olevasta kuvasta näemme myös selvästi tämän autokorrelaation.
R-koodi:

model1<- lm(y[2:100]~x[2:100])
summary(f)
library(lmtest)
bgtest(f,5) 

Model1:n autokorrelaatiokuvat.

Osittaisautokorrelaation piikki ensimmäisessä viiveessä on tyypillistä AR(1)-prosessille, eli dynaamisen regressiomallin muodostaminen on hyvä aloittaa mallilla, jossa on x sekä AR(1)-termi.
R-koodi:

library(forecast)
data <- data.frame(y[2:100],x[2:100])
c <- ts(data,start=c(1999,1),frequency=12) # Muunnan datan aikasarja-muotoon.
arimamodel1 <- Arima(c[,1],xreg=c[,2],order=c(1,0,0)) # Dynaaminen regressiomalli, jossa mukana x ja AR(1)-termi.
Box.test(resid(arimamodel1),type="Ljung",lag=5,fitdf=2) # testaa mallin residuaalien autokorrelaatiota.
res <-resid(arimamodel1)
hist(res, breaks=10, density=10, col="blue", main="Dynaamisen regressiomallin residuaalit", xlab="Residuaalit",freq=FALSE)  # Tekee histogrammin residuaaleista.
x<-seq(min(res),max(res),length=40) 
y<-dnorm(x,mean=mean(res),sd=sd(res)) 
lines(x, y, col="black", lwd=2)

Box-Ljung -testin arvo 0.877 tarkoittaa, että nollahypoteesi (H0: ei autokorrelaatiota virhetermeissä) jää voimaan. Yllä olevasta histogrammista näemme, että residuaalit ovat approksimatiivisesti normaalijakautuneita. Mallin kerroin x:lle on 1.000 ja keskivirhe 0.013, eli pääsimme erittäin lähelle x:n oikeaa arvoa. Testasin vielä muita ARIMA-malleja, mutta tulokset eivät juuri parantuneet. Alla on vielä kuva dynaamisen mallin ennustamista arvoista ja oikeista arvoista.
R-koodi:

plot(arimamodel1$x,col="blue",ylab="y")
lines(fitted(arimamodel1),col="red")
legend("topright",c("Mallin ennusteet","Oikeat y:n arvot"),col=c("red","blue"),lty=c(1,1))