Lab 10: Regresjon med tidsrekke-data II
Johannes Mauritzen
Høsten 2018
I denne labben tar vi opp noen flere emner innenfor regresjon med tidsrekke-data. Her skal vi se på “highly persistent” tidsrekke-data, og i tillegg se litt nærmere på hvordan å modellere tidsrekker som er korrelert gjennom tid. For eksempel, om boligpriser har økt 15 prosent i Oslo i år så er det sansynlig at de kommer til å øke mye neste år også (og året etter det). For å modellere slike serier, skal vi bruke modeller som heter AR (autoregressive), og ARMA(autoregressive Moving Average). I tillegg skal vi se på tidsrekker der variansen endrer seg over tid (ARCH - Autoregressive Conditional Heterodskedascitiy)
Teori
Seriekorrelasjon og Korrelelogram Autokorrelasjonsfunksjon (ACF)
Vi kan bruke ACF til å visualisere hvor mye korrelasjon det er mellom tidspunkter (lags) i vår tidsrekke.
Vi definerer en korrelasjon som følgende:
\(\rho_k = \frac{\gamma_k}{\gamma_0}\)
Der \(\gamma_k\) er kovarians i lag k, og \(\gamma_0\) er varians.
En korrelolgram for serien \(y_t\) vil da vise alle korrelasjonene fra 1 til k:
\(\rho_1 = \frac{cov(y_i, y_{i+1})}{var(y_i)}\)
\(\rho_2 = \frac{cov(y_i, y_{i+2})}{var(y_i)}\)
\(\rho_3 = \frac{cov(y_i, y_{i+3})}{var(y_i)}\)
Vi skaper en seriekorrelert serie (\(y_t = .5*y_{t-1} + u_t\)) og viser en ACF
u=rnorm(100,0,1)
y = (1)
for(t in 2:100){
y[t] = .5*y[t-1] + u[t]
}
t=seq(1:100)
plot(t,y, type="line")## Warning in plot.xy(xy, type, ...): plot type 'line' will be truncated to
## first character
acf(y)
Vi kan tolke dette som ACF plotten som hvor fort en “sjokk” dør ut i en prosess. Her ser det ut som at det tar 5-6 perioder, før en sjokk har sluttet å påvirke en serie.
Vi kan også bruke “partial autocorrelation function.” Disse vil vise autokorrelasjoner mellom ulike lagg, gitt korrelasjonene mellom-laggene.
Vi simulerer \(y_t = .5*y_{t-1} + .3y_{t-2} + u_t\)
u=rnorm(100,0,1)
y2 = c(1, .5*u[1]+u[2])
for(t in 3:100){
y2[t] = .5*y[t-1] + .3*y[t-2] + u[t]
}
t=seq(1:100)
plot(t,y2, type="line")## Warning in plot.xy(xy, type, ...): plot type 'line' will be truncated to
## first character
acf(y2)
pacf(y2)
Seriekorrelasjon
La oss si at vi kjører følgende regresjon:
\(y_t = \alpha + \beta_1 x_t + u_t\)
Hvis feilleddene er autokorrelert (eller seriekorrelert) med hverandre, så kan vi ikke lenger stole på standard feilene fra vår regresjon.
Hvordan kan vi teste for seriekorrelasjon?
For å teste AR(1) seriekorrelasjon, kunne vi ta våre residualer, \(\hat{u_t}\), og kjøre følgende regresjon:
\(\hat{u_t} = \rho \hat{u_{t-1}} + e_t\)
Vi skal se om \(\hat{\rho}\) er signifikant annerledes enn null, og derfor at vi har autokorrelasjon. Men hva er problemet?
Hvis det egentlig er seriekorrelasjon i \(u_t\), så kan vi ikke stole på standard feilene.
Isteden kan vi kjøre
\(\hat{u_t} = \gamma_0 + \gamma_1 x_t + \rho \hat{u_{t-1}}\)
Og se om \(\hat{\rho}\) er annerledes en 0.
For modeller med flere lagger, så kan vi generalisere dette til det som kalles Breusch-Godfrey testen (BG):
Vi kjører
\(\hat{u_t} = \gamma_0 + \gamma_1 x_t + \rho_1 \hat{u_{t-1}} + \rho_2 \hat{u_{t-2}} + ... + e_t\)
Deretter kjører vi en LM-test (som ligner på en F-test) for å se om alle \(\rho\) parameterene sammen er mer er signifikant:
\(LM = (n-q)R^2_{u}\)
der \(n\) er antall datapunkter, \(q\) er antall lag og \(R^2_{u}\) er \(R^2\) statistikken for en “unrestricted” modell, det vil si en som includerer alle laggene.
LM statistikken er da fordelt som en “Chi-square” fordeling.
\(LM \sim X^2_q\)
Data med seriekorrelasjon
Hvis vi finner seriekorrelasjon så har vi flere valgs:
1.) Det beste man kan gjøre er å inkludere dynamikken eksplisit i modellen: sett inn distribuerte lagg og
2.) Justere standard feilene så at de tar hensyn til seriekorrelasjon: dette heter “Newey West” (eller “robust”, eller “HAC”) standard feil.
3.) Noen ganger kan vi transformere vår data ved å differensiere det (se på endringer fra periode til periode). Dette kan også hjelpe å fjerne serie-korrelasjon fra vår data.
library(ggplot2)
chiQ_trekk =data.frame("trekk"= rchisq(100000,2))
perc95 = quantile(chiQ_trekk$trekk, .95)
ggplot(chiQ_trekk, aes(x=trekk)) +
geom_histogram(bins=100) +
geom_vline(xintercept=perc95)
Praksis med R
Noen pakker som vi vil bruke
library(foreign) #for å laste inn data i stata format (read.dta)
library(dynlm) #for dynamiske regresjonsmodeller## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericlibrary(stargazer) #for fine tabeller##
## Please cite as:## Hlavac, Marek (2018). stargazer: Well-Formatted Regression and Summary Statistics Tables.## R package version 5.2.2. https://CRAN.R-project.org/package=stargazerlibrary(ggplot2) #for elegante plotter
library(zoo)
library(reshape2) #for å omformattere data.frame til long-formatt (melt)Vi begynner med å laste inn data om gjennomsnitts-priser og -avkastning på New York Stock Exchange:
nyse <- read.dta("http://fmwww.bc.edu/ec-p/data/wooldridge/nyse.dta")
head(nyse)## price return return_1 t price_1 price_2 cprice cprice_1
## 1 49.75 NA NA 1 NA NA NA NA
## 2 51.44 3.3969820 NA 2 49.75 NA 1.6899990 NA
## 3 52.05 1.1858490 3.3969820 3 51.44 49.75 0.6100006 1.6899990
## 4 52.28 0.4418819 1.1858490 4 52.05 51.44 0.2299995 0.6100006
## 5 54.24 3.7490489 0.4418819 5 52.28 52.05 1.9600030 0.2299995
## 6 53.76 -0.8849619 3.7490489 6 54.24 52.28 -0.4800034 1.9600030Vi formatterer datasette som zoo
nyse_zoo <- zoo(nyse, order.by = nyse$t)
head(nyse_zoo)## price return return_1 t price_1 price_2 cprice cprice_1
## 1 49.75 NA NA 1 NA NA NA NA
## 2 51.44 3.3969820 NA 2 49.75 NA 1.6899990 NA
## 3 52.05 1.1858490 3.3969820 3 51.44 49.75 0.6100006 1.6899990
## 4 52.28 0.4418819 1.1858490 4 52.05 51.44 0.2299995 0.6100006
## 5 54.24 3.7490489 0.4418819 5 52.28 52.05 1.9600030 0.2299995
## 6 53.76 -0.8849619 3.7490489 6 54.24 52.28 -0.4800034 1.9600030og kjører noen dynamiske regresjoner:
return_mod1 <- dynlm(return ~ L(return), data=nyse_zoo)
return_mod2 <- dynlm(return ~ L(return) + L(return, 2), data= nyse_zoo)
return_mod3 <- dynlm(return ~L(return) + L(return, 2) + L(return, 3), data=nyse_zoo)
stargazer(return_mod1, return_mod2, return_mod3, type="text")##
## ===============================================================================
## Dependent variable:
## -----------------------------------------------------------
## return
## (1) (2) (3)
## -------------------------------------------------------------------------------
## L(return) 0.059 0.060 0.061
## (0.038) (0.038) (0.038)
##
## L(return, 2) -0.038 -0.040
## (0.038) (0.038)
##
## L(return, 3) 0.031
## (0.038)
##
## Constant 0.180** 0.186** 0.179**
## (0.081) (0.081) (0.082)
##
## -------------------------------------------------------------------------------
## Observations 689 688 687
## R2 0.003 0.005 0.006
## Adjusted R2 0.002 0.002 0.001
## Residual Std. Error 2.110 (df = 687) 2.112 (df = 685) 2.114 (df = 683)
## F Statistic 2.399 (df = 1; 687) 1.659 (df = 2; 685) 1.322 (df = 3; 683)
## ===============================================================================
## Note: *p<0.1; **p<0.05; ***p<0.01Dette er noen enkle autoregressive regresjoner. “L” står for lag. Det man tester med denne regresjonen er om avkastning på New York Stock Exchange en dag kan si noe om avkastning neste dag og dagen etter det. Dette er definisjonen av persistent. Hvis en serie ikke er persistent, så er det ingen sammenheng - avkastningen en dag er uavhengig av avkastningen dagen før og dagen før det.
Her er det ingen lagget variabler som er signifikant (hvordan vet man det?), og vi kan dermed ikke si at serien er persistent.
Nå laster vi inn pris-data på GE-aksjen
ge_stock <- read.csv("https://jmaurit.github.io/anvendt_macro/data/GE.csv")
head(ge_stock)## Date Open High Low Close Adj.Close Volume
## 1 2012-09-17 21.93 22.05 21.90 22.05 18.52329 79283200
## 2 2012-09-18 21.99 22.24 21.96 22.24 18.68290 38831200
## 3 2012-09-19 22.30 22.49 22.25 22.43 18.84251 40462600
## 4 2012-09-20 22.19 22.47 22.12 22.43 18.98641 43677000
## 5 2012-09-21 22.55 22.69 22.46 22.53 19.07106 66551700
## 6 2012-09-24 22.40 22.45 22.30 22.36 18.92716 36799800ge_stock["Date"]<- as.Date(ge_stock$Date, format="%Y-%m-%d")
ggplot(ge_stock, aes(Date, Close)) +
geom_line()
Vi kan se at tidsrekken av prisene er “highly persistent” - det vil si at prisen i dag er en god indikasjon på hva prisen imorgen kommer til å være i morgen. Tidsrekken er preget av høy auto-korrelasjon og variansen endrer seg over tid. Dette gjør det sånt at man ikke kan bruke OLS på en konsistent måte. Men hvis vi transformerer datarekken til daglig avkastning - det vil si prosent endring per dag, så får vi en datarekke som vi kan bruke.
ge_stock_zoo <- zoo(ge_stock, order.by=ge_stock$Date)
ge_ret <- diff(log(as.numeric(ge_stock_zoo$Close)))
ge_ret_zoo <- zoo(ge_ret, order.by=ge_stock$Date)
plot(ge_ret_zoo)
Da kan vi kjøre noen dynamiske regresjoner. Her er flere AR modeller.
ret_model1 <- dynlm(ge_ret_zoo ~ L(ge_ret_zoo))
ret_model2 <- dynlm(ge_ret_zoo~L(ge_ret_zoo) + L(ge_ret_zoo,2))
ret_model3 <- dynlm(ge_ret_zoo~L(ge_ret_zoo) + L(ge_ret_zoo,2) + L(ge_ret_zoo,3))
stargazer(ret_model1, ret_model2, ret_model3, type="text")##
## =====================================================================================
## Dependent variable:
## -----------------------------------------------------------------
## ge_ret_zoo
## (1) (2) (3)
## -------------------------------------------------------------------------------------
## L(ge_ret_zoo) 0.041 0.043 0.042
## (0.028) (0.028) (0.028)
##
## L(ge_ret_zoo, 2) -0.057** -0.056**
## (0.028) (0.028)
##
## L(ge_ret_zoo, 3) -0.013
## (0.028)
##
## Constant 0.0001 0.0001 0.0001
## (0.0003) (0.0003) (0.0003)
##
## -------------------------------------------------------------------------------------
## Observations 1,257 1,256 1,255
## R2 0.002 0.005 0.005
## Adjusted R2 0.001 0.003 0.003
## Residual Std. Error 0.011 (df = 1255) 0.011 (df = 1253) 0.011 (df = 1251)
## F Statistic 2.095 (df = 1; 1255) 3.057** (df = 2; 1253) 2.102* (df = 3; 1251)
## =====================================================================================
## Note: *p<0.1; **p<0.05; ***p<0.01Her er L en gyldig forkortelse for Lag
Her har vi kjørt tre regresjoner, en AR(1), AR(2) og AR(3) model av avkastninger på GE aksjer. Det vi finner er at koeffisienten på den første laggen er positiv, men ikke signifikant. Den andre laggen er negative og signifikant på 0,05 prosent nivå. F-testen forteller oss at lag1 og lag2 er signifikant sammen. Vi kan tolke koeffisientene på den måten at hvis avkastningen er +100 en dag, så vil det predikere at det blir +4,3 neste dag, men følgende dag går det tilbake -5,7. Serien er kanskje litt persistent: en høy avkastning en dag, kan bety noe høyere avkastning følgende dag, men det blir også en form av “regression to the mean”, at serien deretter går litt ned igjen.
Random Walk
Den enkleste form av en persistent tidsrekke er en random walk. Dette kan skrives:
\[y_t = y_{t-1} + e_t\]
Intuitivt - prisen i dag er lik prisen i går pluss en tilfeldig (“stokastisk”) hopp, som vi representerer med \(e_t\). Random walk og noen litt mer kompliserte modeller blir brukt til å modellere oljepriser, aksjepriser, og også noen makroøkonomiske variabler.
Vi kan simulere random walk enkelt med r:
#lager en matrise som er 50x30
y<-matrix(,nrow=50, ncol=30)
#fyller den ved bruke av en for-loop
for(s in 1:30){
e <- rnorm(50) #her er vår stokastisk variabel - vi simulerer 50 trekk fra en normal-fordeling
y[,s]<-cumsum(e) #y er da modellert som summen av alle "sjokkene", e.
}
y = data.frame(y)
y["t"]<-1:50
y_l=melt(data.frame(y), id.var="t")
ggplot(y_l, aes(x=t, y=value, group=variable)) +
geom_line(alpha=.5)
I figuren vises 30 simulerte tidsrekker av 50 perioder. Vi kan se at variansen øker over tid og at serien ikke beveger seg mot en bestemt gjennomsnittsverdi - noe som gjør OLS ugyldig. Serien er det vi kaller for ikke-stasjonær (som vi kommer tilbake til).
I beregningen, så regner vi ut \(y_t\) som en kumulativt sum av \(e_t\) ’ene. det vil si \(y_t = e_0 + e_1 + ... + e_{t-1} + e_t\). Kan du forklare hvorfor?
En annen form av random walk er random walk med drift. Det vil si at random walk mønsteret har en bestemt retning. Dette kan vi skrive som
\[y_t = \alpha_0 + y_{t-1} + e_t\]
#lager en matrise som er 50x30
y<-matrix(,nrow=50, ncol=30)
#fyller det ved bruke av en for-loop
for(s in 1:30){
e <- rnorm(50) #her er vår stokastisk variabel - vi simulerer 50 trekk fra en normal-fordeling
y[,s]<-cumsum(2+e) #y er da modellert som summen av 2 og alle "sjokkene", e.
}
y = data.frame(y)
y["t"]<-1:50
y_l=melt(data.frame(y), id.var="t")
ggplot(y_l, aes(x=t, y=value, group=variable)) +
geom_line(alpha=.5)
Man kunne bruke random walk med drift til å modellere aksjepriser over tid, for eksempel. Vi ser at serien har en tydelig trend oppover, og at den har den vanlige egenskapen til en en random walk - variansen som øker over tid og en ikke-definert gjennomsnitt (forventning).
Differanser av highly persistent data
Highly persistent data (random walk, osv) er, som nevnt, vanskelig å jobbe med i regresjons-sammenheng. Heldigvis er det ofte en enkel løsning - differensiering (first difference).
Tidsrekken:
\[y_t = \alpha_0 + y_{t-1} + e_t\]
blir:
\[\Delta y_t = y_t- y_{t-1} = \alpha + e_t\]
Vi kan simulere:
for(s in 1:30){
e <- rnorm(50) #her er vår stokastisk variabel - vi simulerer 50 trekk fra en normal-fordeling
y[,s]<-cumsum(2+e) #y er da modellert som summen av 2 og alle "sjokkene", e.
}
y = data.frame(y)
y["t"]<-1:50
y <- zoo(y, order.by = y$t)
Dy <- diff(y)
Dy <- data.frame(Dy)
Dy["t"]<-1:49
Dy_l=melt(Dy, id.var="t")
ggplot(Dy_l, aes(x=t, y=value, group=variable)) +
geom_line(alpha=.2) +
geom_hline(yintercept=2)
Nå kan vi se at variansen av datasettet ser omtrent konstant gjennom tid. Vi kunne klart å beregne en OLS regresjon som ville gitt oss omtrent en vannrett linje.
Regresjon med first-difference
Vi bruker fruktbarhetsdataen igjen for å vise hvordan å bruke first-difference:
fertil3 <- read.dta("http://fmwww.bc.edu/ec-p/data/wooldridge/fertil3.dta")
fertil3_zoo <- zoo(fertil3, order.by=fertil3$year)
plot(fertil3_zoo$gfr)
plot(diff(fertil3_zoo$gfr))
Datarekken for fruktbarhet (gfr) ser ut som den er persistent, men vi kan ta en first difference.
Da kan vi kjøre regresjon av fruktbarhet på skattefradrag (pe).
fert_d_mod1 <- dynlm(diff(gfr) ~ diff(pe), data = fertil3_zoo)
fert_d_mod2 <- dynlm(diff(gfr) ~ diff(pe) + L(diff(pe)) + L(d(pe), 2) , data = fertil3_zoo)
stargazer(fert_d_mod1, fert_d_mod2, type="text")##
## ============================================================
## Dependent variable:
## ----------------------------------------
## diff(gfr)
## (1) (2)
## ------------------------------------------------------------
## diff(pe) -0.043 -0.036
## (0.028) (0.027)
##
## L(diff(pe)) -0.014
## (0.028)
##
## L(d(pe), 2) 0.110***
## (0.027)
##
## Constant -0.785 -0.964**
## (0.502) (0.468)
##
## ------------------------------------------------------------
## Observations 71 69
## R2 0.032 0.232
## Adjusted R2 0.018 0.197
## Residual Std. Error 4.221 (df = 69) 3.859 (df = 65)
## F Statistic 2.263 (df = 1; 69) 6.563*** (df = 3; 65)
## ============================================================
## Note: *p<0.1; **p<0.05; ***p<0.01Oppgave: Hvordan kan man tolke resultatene? Hva blir resultatene om man ikke bruker first-diff?
Svar:
Når vi gjør en first-difference, så spør vi – Hva er effekten av en endring i skattefradrag på endring i fruktbarhet? Eller – Hvis vi endrer skattefradrag, hva blir endringen i fruktbarhet? Koeffisienten på skattefradrag og laggen er negativ, men koeffisienten er ikke signifikant annerledes enn 0. Vi kan, med andre ord, ikke konkludere med noe.
Hvis vi ikke brukte first-difference, ville vi gjort som følgende:
fert_mod2 <- dynlm(gfr ~ pe + L(pe) + L(pe, 2) , data = fertil3_zoo)
summary(fert_mod2)##
## Time series regression with "zoo" data:
## Start = 1915, End = 1984
##
## Call:
## dynlm(formula = gfr ~ pe + L(pe) + L(pe, 2), data = fertil3_zoo)
##
## Residuals:
## Min 1Q Median 3Q Max
## -30.146 -16.062 -0.127 17.562 31.842
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 93.15791 4.49965 20.703 <2e-16 ***
## pe -0.01584 0.14026 -0.113 0.910
## L(pe) -0.02134 0.21523 -0.099 0.921
## L(pe, 2) 0.05390 0.13811 0.390 0.698
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 19.78 on 66 degrees of freedom
## Multiple R-squared: 0.006135, Adjusted R-squared: -0.03904
## F-statistic: 0.1358 on 3 and 66 DF, p-value: 0.9383I denne tilfellen er det ikke store kvalitative forskjeller mellom de to regresjonene – vi får negative koeffisienter, men disse er ikke signifikante annerledes enn 0. Dette kan tyde på at vi ikke behvøde å differensiere de opprinnelige seriene for å få valide inferens.
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 standard feilene er ikke lenger “efficient” - det vil si at vi ikke kan stole på standard feilene og p-verdier. Det er best å prøve å kontrollere for seriekorrelasjon.
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-05Her 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. (Man kan kjøre Breusch-godfrey testen ved å bruke bgtest funksjonen i lmtest pakken)
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.02271Her 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.7752Ved 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-05Vi 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)
library(lmtest)
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.002264En 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 0dwtest(phillips_mod2)##
## Durbin-Watson test
##
## data: phillips_mod2
## DW = 1.7696, p-value = 0.1783
## alternative hypothesis: true autocorrelation is greater than 0Sist 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”) (Se Hill, et al. k. 16).
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.76Litt 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.