library(foreign)
library(zoo)
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
library(dynlm)

Seriekorrelasjon

Vi har sett på seriekorrelasjon/autokorrelasjon før - at feilleddet er korrelert over tid i en regresjon. Hvis alt annet er bra (konstant varians, osv) så påvirker ikke dette om estimatet er biased, men standardavvik er ikke lenger “efficient” - det vil si at vi ikke kan stole på standard avvik og p-verdier. Det er best å prøve å kontrollere for seriekorrelasjon.

(ps Det finnes også noen videor på nasjonalkursrom som dekker seriekorrelasjon, hvis du ønsker en gjennomgang fra et litt annet perspektiv.)

Testing for seriekorrelasjon - Phillips curve

Phillipskurven er den observerte forholdet mellom inflasjon og arbeidsledighet. Vi tar inn data for å prøve å måle dette forholdet

phillips = read.dta("http://fmwww.bc.edu/ec-p/data/wooldridge/phillips.dta")
phillips_zoo = zoo(phillips, order.by=phillips$year)

#scatter plot
plot(phillips[c("unem", "inf")])

#line plots
plot(phillips_zoo$unem)

plot(phillips_zoo$inf)

Vi kjører en enkel regresjon og da får vi ut residualene (husker du hvordan vi definerer residualer?). Vi skal plotte residualene i mot “fitted” verdiene (det vil si, prediksjonene som modellen gir for inflasjon gitt arbeidslediget.) Hvis modellen er grei, burde vi ikke se noe mønster her. Hvis vi ser en tydelig mønster, så er det noe vi burde prøve å modellere. Dette inkluderer seriekorrelasjon – ser det ut som at en høy verdi er følgt av en annen høy verdi, osv?

phillips_mod = dynlm(inf ~ unem, data=phillips_zoo)
phill_resid = resid(phillips_mod)
phill_fitted = fitted(phillips_mod)

#fitted vs. residual plott
plot(phill_fitted, phill_resid)

Fitted vs. residual plotten ser egentlig ganske bra ut (det burde se ganske tilfeldig - ingen bestemte trender.)

Men vi kan teste for seriekorrelasjon (autokorrelasjon) også på en mer formelt måte:

ar_test = dynlm(phill_resid ~ L(phill_resid))
summary(ar_test)
## 
## Time series regression with "zoo" data:
## Start = 1949, End = 1996
## 
## Call:
## dynlm(formula = phill_resid ~ L(phill_resid))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -8.0764 -1.2577 -0.2327  0.9657  6.5375 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     -0.1134     0.3594  -0.316    0.754    
## L(phill_resid)   0.5730     0.1161   4.934  1.1e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.49 on 46 degrees of freedom
## Multiple R-squared:  0.346,  Adjusted R-squared:  0.3318 
## F-statistic: 24.34 on 1 and 46 DF,  p-value: 1.098e-05

Her kjører vi en regresjon av residualene på sin lagg. Siden laggen er signifikant, har vi grunn til å tro at residualene er seriekorrelert. Dette er en enkel form av det som er kalt Breusch-Godfrey testen

Vi kan gjøre det samme med en model som er “expectations augmented”, det vil si at forholdet er mellom endring i inflasjon fra år til år og arbeidsledighet.

Oppgave:

Kjør en regresjon av en “expectations augmented” phillips-kurve, tolk resultatene. Kan du knytte det til vår diskusjon om phillips-kurven tidligere i semesteret? Test om det finnes seriekorrelasjon i feilleddet fra denne modellen.

Svar

phillips_mod2 = dynlm(d(inf) ~ unem, data=phillips_zoo)
summary(phillips_mod2)
## 
## Time series regression with "zoo" data:
## Start = 1949, End = 1996
## 
## Call:
## dynlm(formula = d(inf) ~ unem, data = phillips_zoo)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -9.1293 -0.7820  0.0611  0.9887  5.3600 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)   3.0306     1.3768   2.201   0.0328 *
## unem         -0.5426     0.2302  -2.357   0.0227 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.451 on 46 degrees of freedom
## Multiple R-squared:  0.1078, Adjusted R-squared:  0.0884 
## F-statistic: 5.558 on 1 and 46 DF,  p-value: 0.02271

Her ser vi at endringer i inflasjon er signifikant og negativt relatert til arbeidsledighet, som teorien om den “expectations-augmented” Phillipskurven ville forutsi.

Vi kan videre teste for seriekorrelasjon:

phill_resid2 = resid(phillips_mod2)
ar_test2 = dynlm(phill_resid2 ~ L(phill_resid2))
summary(ar_test2)
## 
## Time series regression with "zoo" data:
## Start = 1950, End = 1996
## 
## Call:
## dynlm(formula = phill_resid2 ~ L(phill_resid2))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7.4062 -0.8323 -0.0874  0.7605  5.2493 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)
## (Intercept)      0.19417    0.30038   0.646    0.521
## L(phill_resid2) -0.03559    0.12389  -0.287    0.775
## 
## Residual standard error: 2.059 on 45 degrees of freedom
## Multiple R-squared:  0.001831,   Adjusted R-squared:  -0.02035 
## F-statistic: 0.08254 on 1 and 45 DF,  p-value: 0.7752

Ved bruk av pakken lmtest kan vi teste for autokorrelasjon (AR(1)) med Breusch-Godfrey-testen med bare en linje med kode.

library(lmtest)
bgtest(inf ~ unem, data=phillips_zoo)
## 
##  Breusch-Godfrey test for serial correlation of order up to 1
## 
## data:  inf ~ unem
## LM test = 18.472, df = 1, p-value = 1.725e-05

Vi kan også teste for seriekorrelasjon med lengre virkning

\[e_t = \alpha_1 e_{t-1} + \alpha_2 e_{t-2} + ... + \alpha_n e_{n} \]

Vi bruker barium datasettet til å se på dette:

barium =  read.dta("http://fmwww.bc.edu/ec-p/data/wooldridge/barium.dta")
barium_zoo = zoo(barium, order.by = barium$t, frequency=12)

barium_mod = dynlm(log(chnimp) ~ log(chempi) + log(gas) + log(rtwex) + befile6 + affile6 + afdec6, data=barium_zoo)
bgtest(barium_mod, order=3, type="F")
## 
##  Breusch-Godfrey test for serial correlation of order up to 3
## 
## data:  barium_mod
## LM test = 5.1247, df1 = 3, df2 = 121, p-value = 0.002264

En annen test som brukes ofte for seriekorrelasjon er durbin-watson. Vi viser det med phillips-kurv datasettet:

#funksjonen dwtest kommer med pakken lmtest
dwtest(phillips_mod) 
## 
##  Durbin-Watson test
## 
## data:  phillips_mod
## DW = 0.8027, p-value = 7.552e-07
## alternative hypothesis: true autocorrelation is greater than 0
dwtest(phillips_mod2)
## 
##  Durbin-Watson test
## 
## data:  phillips_mod2
## DW = 1.7696, p-value = 0.1783
## alternative hypothesis: true autocorrelation is greater than 0

Sist men ikke minst, kan vi se på autokorrelasjonsfunksjoner (ACF) og “partial” autokorrelasjonsfunksjoner (pacf)

Vi kan se på inflasjon serien i vår phillipskurv datasett og bruke en autokorrelasjonsfunksjon (som visuelt heter “Correlelogram”).

acf(phillips_zoo$inf)

Du kan tenke på en ACF som om hvis det er en type “sjokk” i t=0 av 1, hvor lenge vedvarer denne sjokken. Kanskje vi får en liten sjokk i inflasjon – høyere bensinpriser på grunn av orkanen i USA – hvor lenge påvirker det inflasjonsserien? Vi ser at inflasjon har noe seriekorrelasjon, men at en sjokk dør ut etter 4-5 år.

Vi kan også se på en pACF:

pacf(phillips_zoo$inf)

Vi kan tolke en pACF mer direkt som hvor mange lagger ville være signifikant i en regresjon. Her ser det ut som det er hovedsakelig en lag som signifikant, og at vi kunne modellere inflasjon som en AR(1) modell.

Regresjon som er robust til seriekorrelasjon

Vi begynner med å laste inn data om Brent oljepriser fra FRED data-basen:

brent=read.csv("https://fred.stlouisfed.org/graph/fredgraph.csv?chart_type=line&recession_bars=on&log_scales=&bgcolor=%23e1e9f0&graph_bgcolor=%23ffffff&fo=Open+Sans&ts=12&tts=12&txtcolor=%23444444&show_legend=yes&show_axis_titles=yes&drp=0&cosd=1987-05-01&coed=2016-08-01&height=450&stacking=&range=&mode=fred&id=MCOILBRENTEU&transformation=lin&nd=1987-05-01&ost=-99999&oet=99999&lsv=&lev=&mma=0&fml=a&fgst=lin&fgsnd=2009-06-01&fq=Monthly&fam=avg&vintage_date=&revision_date=&line_color=%234572a7&line_style=solid&lw=2&scale=left&mark_type=none&mw=2&width=1168")
head(brent)
##         DATE MCOILBRENTEU
## 1 1987-05-01        18.58
## 2 1987-06-01        18.86
## 3 1987-07-01        19.86
## 4 1987-08-01        18.98
## 5 1987-09-01        18.31
## 6 1987-10-01        18.76

Litt rengjøring…

brent["DATE"] = as.Date(brent$DATE)
brent_ts = ts(brent$MCOILBRENTEU, start=c(1987,1), end=c(2016,8), frequency=12)
plot(brent_ts)

Oppgave

  • Hva kan man si om denne serien bare ved å se på den?

  • Bruk ACF og pACF til å analysere hvor autokorrelert denne serien er. Hva sier dette intuitivt?

  • Hvis det er en høy grad av seriekorrelasjon, kan serien være “ikke-stasjonært”, som gjør det vanseklig å analysere med OLS. En løsning er å ta first-difference, eller å se på avkastning. Gjør dette med oljeprisdata og kjør en ACF igjen. Hva har skjedd?

  • Hvis det er fortsatt autokorrelasjon i serien, kjør en dynamisk regresjon som modellerer denne seriekorrelasjonen eksplisitt.

### Svar
#En ting vi kan se er variansen virker som at det øker over tid -- større svingninger i oljeprisen etter hvert. Vi kan også se at serien ikke går tilbake til en viss gjennomsnittsverdi. Vi må dermed være forsiktig med å kjøre OLS. 


acf(brent_ts)

#Om det ikke var significant autokorrelasjon, ville de vertikale linjene vært nær 0. Isteden ser vi at det er mye autokorrelasjon i datasettet. Fra det vi så tidligere, så kan oljeprisen minne på en random walk - den er highly persistent. 

#La oss transformere datasette til månedlig avkastning og sjekke acf igjen. 

brent_ret = diff(log(brent_ts))
plot(brent_ret)

acf(brent_ret)

#Nå ser det bedre ut, men det ser ut som avkastning har fortsatt autokorrelasjon. Vi kan modellere dette med *dynlm* pakken


brent_mod = dynlm(brent_ret~ L(brent_ret) + L(brent_ret, 2))
summary(brent_mod)
## 
## Time series regression with "ts" data:
## Start = 1987(4), End = 2016(8)
## 
## Call:
## dynlm(formula = brent_ret ~ L(brent_ret) + L(brent_ret, 2))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.90760 -0.05116  0.00785  0.05934  0.43185 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     -0.0001582  0.0053258  -0.030 0.976313    
## L(brent_ret)     0.2065738  0.0534559   3.864 0.000133 ***
## L(brent_ret, 2) -0.0090583  0.0534737  -0.169 0.865581    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1001 on 350 degrees of freedom
## Multiple R-squared:  0.04199,    Adjusted R-squared:  0.03652 
## F-statistic: 7.671 on 2 and 350 DF,  p-value: 0.000549
#Den første laggen er positiv og signifikant. Vi kan dermed tolke det til å bety at en endring i oljeprisen en periode vil ofte føre til, i gjennomsnitt, en økning i neste periode også.

Robust t-test

Når vi finner seriekorrelasjon, så er den beste løsningen å prøve å modellere autokorrelasjonen direkt i modellen, som vi gjorde i eksempelen over. Men det er også mulig å estimere en regresjon der t-testene er robust til seriekorrelasjon.

I r kan vi gjøre det ved å ta i bruk pakken sandwich og ta i bruk det som heter en “Newey West” estimering, som tar hensyn til seriekorrelasjon i leddet.

library(sandwich)
brent_mod2 = dynlm(brent_ret~ L(brent_ret))

coeftest(brent_mod2, vcov =  NeweyWest)
## 
## t test of coefficients:
## 
##                 Estimate  Std. Error t value Pr(>|t|)  
## (Intercept)  -2.0647e-05  5.1383e-03 -0.0040  0.99680  
## L(brent_ret)  2.0491e-01  8.7745e-02  2.3353  0.02009 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

I kommandoen som vi ser over så representerer “HC” Heteroskedacity and SerialCorrelation Robust estimatore

ARCH

ARCH står for Autoregressive Conditional Heteroscedasticity. Denne type modell tar hensyn til en viss type heteroskedasitet (ikke-konstant varianse) i dataen. Man finner ofte disse modellene anvendt i finans.

Hvis vi har en modell:

\[y_t = \beta_0 + \beta_1 x_{t1} + \beta_2 x_{t2} + ... + \beta_k x_{tk}+ u_t\]

Feilleddet følger en ARCH prossess hvis:

\[E(u^2_t|u_{t-1}, u_{t-2},...) = \alpha_0 + \alpha_1 u^2_{t-1}\]

Variansen i t er avhengig av variansen i t-1, og vi kan modellere dette via OLS.

Som en eksempel, ser vi på avkastning på aksjer:

nyse = read.dta("http://fmwww.bc.edu/ec-p/data/wooldridge/nyse.dta")
head(nyse)

Vi omformatterer til tidsrekke-data via zoo og kaster vekk observasjonene der man mangler data (NA).

nyse_zoo = zoo(nyse, order.by = nyse$t)

#få vekk NA data
nyse_zoo = nyse_zoo[!is.na(nyse_zoo$return), ]

Vi begynner med å kjøre en vanlig AR(1) regresjon av avkastningningen. Vi tar deretter kvadraten av residualene fra denne regresjonen og kjører en AR(1) regresjon på dem.

arch_mod1 = dynlm(return ~ L(return), data=nyse_zoo)

#kvadrat-residual
residual.sq = resid(arch_mod1)^2

#modell for kvadrat-residual
ARCHreg = dynlm(residual.sq ~ L(residual.sq)-1)

coeftest(ARCHreg)
## 
## t test of coefficients:
## 
##                Estimate Std. Error t value  Pr(>|t|)    
## L(residual.sq) 0.424453   0.034545  12.287 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Her ser vi at laggen av kvadrat-residualen er signifikant. Det vil si at variansen i denne modellen er avhengig av tidligere varians. Noen viktige ting vi kan si om dette:

  • En av våre forutsetninger for å bruke OLS med tidsrekkedata var at variansen er konstant. Her har vi data der dette er ikke tilfelle. Variansen varierer over tid. For å få riktig inferens må vi derfor ta hensyn til dette.

  • At en serie følger en ARCH-prossess er også interessant i seg selv. For denne serien sier det at perioder med høy varians (høy risiko) er persistent. Hvis det var høy varians i en periode, så er det sansynlig at det blir høy varians neste periode også. Dette kan bety, for eksempel, at når det først er frykt på markedet, så kan det vedvare og føre til mer frykt.