I den forrige labben, gikk vi gjennom hva stasjonæritet er, hvorfor det er viktig at vi har stasjonære serier og vi så på noen strategier for å transformere ikke-stasjonære serier. I denne labben, fortsetter vi vår diskusjon av stasjonæritet. Vi skal se på hvordan å formelt teste for ikke-stasjonære serier og se på noen eksempler av å kjøre regresjon med ikke-stasjonær data.

Teste for I(1)

Hvordan vet vi at vår serie er ikke-stasjonær? Fram til nå, har vi sett på våre serier til å vurdere om de er stasjonær eller ikke. Men det ville være fint å ha en litt mer treff-sikker måte å etablere om serien er stasjonær eller ikke.

Hvis vi ser på en serie som vi tror er AR(1), så kunne vi estimere regresjonen:

\(y_t = \rho y_{t-1} + u_t\)

Vi kan simulere data og prøve:

Jeg skal simulere en tidsrekke fra en AR(1) modell der jeg velger \(\rho=0,8\). Deretter skal jeg kjøre en regresjon av \(y_t\) på sin lag \(y_{t-1}\) for å få et estimat av \(\rho\) som jeg kaller \(\hat{\rho}\).

#først setter vi rho=.8
library(zoo)
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
library(dynlm)

u=rnorm(100,0,1)
y = (1) #første verdi
rho = .8 #AR ledd
for(t in 2:100){
  y[t] = rho*y[t-1] + u[t] 
}
t=seq(1:100)

y_zoo = zoo(y, order.by=t)

mod1 = dynlm(y_zoo ~ lag(y_zoo)-1)
summary(mod1)
## 
## Time series regression with "zoo" data:
## Start = 1, End = 99
## 
## Call:
## dynlm(formula = y_zoo ~ lag(y_zoo) - 1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.63376 -0.74851  0.06194  0.79333  3.03639 
## 
## Coefficients:
##            Estimate Std. Error t value Pr(>|t|)    
## lag(y_zoo)  0.78306    0.06314    12.4   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9538 on 98 degrees of freedom
## Multiple R-squared:  0.6108, Adjusted R-squared:  0.6068 
## F-statistic: 153.8 on 1 and 98 DF,  p-value: < 2.2e-16

Vi burde få en koeffisient på \(y_{t-1}\) (som heter lag(y)) som er ganske nær 0,8 (forhåpentligvis), og ihvertfal noe som indikerer at \(\rho<1\). Men hva skjer om \(\rho=1\)

#nå setter vi rho=1

u=rnorm(100,0,1)
y = (1) #første verdi
rho = 1 #AR ledd
for(t in 2:100){
  y[t] = rho*y[t-1] + u[t] 
}
t=seq(1:100)

y = zoo(y, order.by=t)

mod2 = dynlm(y ~ lag(y))
summary(mod2)
## 
## Time series regression with "zoo" data:
## Start = 1, End = 99
## 
## Call:
## dynlm(formula = y ~ lag(y))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.43402 -0.77889 -0.08087  0.70981  2.19807 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.21348    0.18545   1.151    0.253    
## lag(y)       0.94957    0.03336  28.465   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.044 on 97 degrees of freedom
## Multiple R-squared:  0.8931, Adjusted R-squared:  0.892 
## F-statistic: 810.3 on 1 and 97 DF,  p-value: < 2.2e-16

Hvis man kjørte denne regresjonen flere ganger med ulike simulasjoner av data, det man ville oppleve var at resultatene spriker. Noen ganger får man et estimat som er nær 1, andre ganger får man estimater som er langt mindre eller større.

Problemet med å kjøre en slik regresjon for å teste stasjonæritet er at når serien ikke er stasjonær, så er t-testene ugyldig.

Vi kommer tilbake til hoved-beskjeden om stasjonæritet: når vi har ikke-stasjonær data, så kan vi ikke stole på regresjonsresultatene. Derfor kan vi heller ikke bruke vanlig regresjon til å teste for stasjonæritet.

Vi bruker heller noe som heter en Dickey-Fuller test.

Dickey-Fuller test:

Vi igjen starter med vår enkel AR(1) modell:

\[y_t = \rho y_{t-1} + u_t\]

Null-hypotesen vår er:

\(H_0: \rho=1\) ikke-stasjonær

mot alternativet

\(H_A: \rho<1\)

Vi husker at hvis vi transformerer til endringer, så vil det ofte gjøre en stasjonær serie ikke-stasjonær. La oss gjøre det med vår AR(1) modell:

\(y_t - y_{t-1} = \rho y_{t-1} + e_t - y_{t-1}\)

\(\Delta y_{t-1} = (\rho - 1) y_{t-1} + u_t\)

Dermed, hvis vi definerer \(\theta = \rho-1\), får vi følgende regresjonsmodel:

\(\delta y_t = \theta y_{t-1} + u_t\)

Nå kan vi teste for stasjonæritet ved følgende hypoteser:

\(H_0: \theta=0\): ikke stasjonær (null hypotese)

\(H_A: \theta<0\): stasjonær (alternativ hypotese)

Hvis vi utfører denne testen med en spesiell fordeling som heter Dickey-Fuller fordelingen, så utgjør det Dickey-fuller testen for stasjonæritet.

Augmented Dickey Fuller test

Ofte, må man ta hensyn til mer dynamikk og seriekorrelasjon i \(\Delta y_t\) serien, da bruker vi det vi kaller for en Augmented Dickey Fuller Test (ADF-testen).

En regresjon med AR(1) seriekorrelasjon i \(\Delta y_t\) kan vi skrive som følgende:

\(\Delta y_t = \theta y_{t-1} + \gamma \Delta y_{t-1} + e_t\)

Og generelt med p-lagger:

\(\Delta y_t = \theta y_{t-1} + \gamma_1 \Delta y_{t-1} + ... + \gamma_p \Delta y_{t-p} + e_t\)

og deretter teste om \(\theta = 0\) med Dickey-fuller fordelingen som vanlig.

Case study: oljepris og kronekursen

I denne casen skal vi bruke noen av verktøyene vi har lært for å analysere hvordan oljeprisen utvikler seg - økonometrisk-sett, og hvordan den påvirker noen aspekter av den norske økonomien.

library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.0 ──
## ✓ ggplot2 3.3.3     ✓ purrr   0.3.4
## ✓ tibble  3.0.6     ✓ dplyr   1.0.4
## ✓ tidyr   1.1.2     ✓ stringr 1.4.0
## ✓ readr   1.3.1     ✓ forcats 0.5.0
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
library(foreign)
library(dynlm)
library(zoo)

Figuren viser data om Brent oljepris som er hentet inn fra St. Louis Reserve Bank FRED databasen:

oil_price_series = 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=2018-09-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")
brent = oil_price_series
colnames(brent) = c("date", "brentprice")
brent["date"] = as.Date(brent$date)
ggplot(brent, aes(date, brentprice)) +
  geom_line()

Hvis vi skal vurdere stasjonæritet visuelt, ville du sagt at serien ser stasjonær ut?

La oss først kjøre en naive AR(1) regresjon av vår data. Resultatene kan ihvertfall gi oss en indikasjon på stasjonæritet.

brent_ts = zoo(brent$brentprice, order.by=brent$date)
I1_mod = dynlm(brent_ts~ L(brent_ts))
summary(I1_mod)
## 
## Time series regression with "zoo" data:
## Start = 1987-06-01, End = 2018-09-01
## 
## Call:
## dynlm(formula = brent_ts ~ L(brent_ts))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -25.366  -1.483   0.135   1.696  14.116 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 0.554119   0.391584   1.415    0.158    
## L(brent_ts) 0.991385   0.006941 142.822   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.451 on 374 degrees of freedom
## Multiple R-squared:  0.982,  Adjusted R-squared:  0.9819 
## F-statistic: 2.04e+04 on 1 and 374 DF,  p-value: < 2.2e-16

Vi estimerer en koeffisient ganske nær 1, så det lover dårlig.

Nå kjører vi en ADF test. Du kan gjøre dette automatisk i både Stata og R.

I Stata bruker man kommandoen dfuller

I R bruker man adf.test i pakken tseries:

library(tseries)
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
adf.test(brent_ts)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  brent_ts
## Dickey-Fuller = -2.4318, Lag order = 7, p-value = 0.3948
## alternative hypothesis: stationary

Tolkning:

I vår Augmented Dickey-Fuller testen er alternative hypotesen at serien er stasjonær. Vi får en Dickey-Fuller test-statistikk på -2.43 og p-verdi av 0,3948. Det vil si at under vår null-hypotese av en ikke-stasjonær serie, ville vi fått test-statistikk som er lavere (mer negativ) 39% av tiden. Hvis vi har en cut-off verdi av 5%, så er vi langt i fra å kunne forkaste vår null-hypotese. Vi burde dermed anta at vi har en stasjonær serie.

Men vi kan gjøre serien om til avkastninger (prosent endringer) og se om det fikser problemet.

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

Denne ser mer stasjonær ut - den har en noenlunde konstant varians og en gjennomsnittsverdi som serien beveger seg rundt, men vi kan sjekke formelt med nok en adf-test:

#vi kvitter oss med manglede (NA) data
brent_ret = brent_ret[!is.na(brent_ret)] 

#dickey-fuller test
adf.test(brent_ret, k=1)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  brent_ret
## Dickey-Fuller = -12.071, Lag order = 1, p-value = 0.01
## alternative hypothesis: stationary

Her ser vi at Dickey-Fuller testen er -12,07 med p-verdi av 0,01. Det betyr at hvis null-hypotesen er riktig (ikke-stasjonær), så vil vi få en test-statistikk av mindre (mer negativ) enn -12,07 1% av tiden. Hvis vi har en cut-off verdi av 5%, så kan vi forkaste vår null-hypotese. Vi tror nå at vår tidsrekke er stasjonær.

Modellering av oljepris

Vi kan begynne vår analyse av oljeprisen på å undersøke dynamikken av oljeprisen, i avkastning form.

Vi kjører en AR(3) modell (regresjon av opljepris-avkastningen på 3 lagger.)

ac_mod = dynlm(brent_ret ~ L(brent_ret) + L(brent_ret,2) + L(brent_ret, 3))
summary(ac_mod)
## 
## Time series regression with "zoo" data:
## Start = 1987-09-01, End = 2018-09-01
## 
## Call:
## dynlm(formula = brent_ret ~ L(brent_ret) + L(brent_ret, 2) + 
##     L(brent_ret, 3))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.27302 -0.05347  0.00531  0.05503  0.41838 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      0.0029983  0.0044636   0.672    0.502    
## L(brent_ret)     0.2720204  0.0520948   5.222 2.97e-07 ***
## L(brent_ret, 2) -0.0328973  0.0539450  -0.610    0.542    
## L(brent_ret, 3) -0.0006892  0.0520807  -0.013    0.989    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.08605 on 369 degrees of freedom
## Multiple R-squared:  0.07027,    Adjusted R-squared:  0.06271 
## F-statistic: 9.296 on 3 and 369 DF,  p-value: 6.117e-06

Ja, det ser ut som den første laggen er signifikant.

Vi kunne også bruke ACF og pACF figurer til å undersøke dynamikken i vår serie

Kronekursen og oljeprisen

Nå skal vi prøve å modellere effekten oljepris har på kronekursen.

Vi bruker data om valutakurs fra Norges Bank:

#Here is the cleaning of the data
valutakurs = read_csv2("https://data.norges-bank.no/api/data/EXR/M..NOK.SP?startPeriod=2000-09-15&endPeriod=2020-09-15&format=csv&locale=no")

kronekurs = valutakurs %>% filter(BASE_CUR =="USD")
kronekurs = kronekurs %>% select("TIME_PERIOD", "OBS_VALUE")
colnames(kronekurs) = c("date", "NOK_USD")
dates = as.yearmon(kronekurs$date, format="%Y-%m")
kronekurs["date"] = as.Date(dates)
kronekurs %>% write_csv("data/kroneKurs.csv")
kronekurs = read_csv("https://jmaurit.github.io/anvendt_macro/lab/data/kroneKurs.csv")
## Parsed with column specification:
## cols(
##   date = col_date(format = ""),
##   NOK_USD = col_double()
## )
#kronekurs = read_csv("data/kroneKurs.csv")

#få valutakurs-formatt
kronekurs["USD_NOK"] = 1/kronekurs["NOK_USD"]

#bruker month og year, og funksjonen merge til å sette sammen seriene
oil_krone = merge(brent, kronekurs, id="date")
head(oil_krone)
##         date brentprice NOK_USD   USD_NOK
## 1 2000-10-01      30.96  9.3613 0.1068228
## 2 2000-11-01      32.55  9.3369 0.1071019
## 3 2000-12-01      25.66  9.0662 0.1102998
## 4 2001-01-01      25.62  8.7784 0.1139160
## 5 2001-02-01      27.50  8.9117 0.1122120
## 6 2001-03-01      24.50  8.9742 0.1114305
#Klassifiserer oljeprisen og kronekurs som zoo-serie
oil_krone_zoo = zoo(oil_krone[c("brentprice", "USD_NOK")], order.by=oil_krone$date) 

plot(x=oil_krone_zoo)

Det første vi burde se er at våre serier er stasjonære. Vi har testet dette allerede med vår oljepris data. Vi kan gjøre det samme med kronekursen.

library(tseries)
adf.test(oil_krone_zoo$USD_NOK)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  oil_krone_zoo$USD_NOK
## Dickey-Fuller = -1.97, Lag order = 5, p-value = 0.5886
## alternative hypothesis: stationary

Her får vi en Dickey-Fuller statistikk som er -2.57 og en p-verdi av .33. Det betyr at under vår null-hypotese (der vi antar ikke-stasjonæritet), så er det en 33% sjanse at man får -2.57 eller lavere. Dette er ikke godt nok for oss å forkaste null-hypotesen (ikke-stasjonær). Vi må gå ut i fra at denne serien også er ikke-stasjonær.

Vi transfomerer derfor også vår valutakurs data til avkastninger.

oil_krone["oil_price_ret"] = log(oil_krone$brentprice) - lag(log(oil_krone$brentprice))
oil_krone["USD_NOK_ret"] = log(oil_krone$USD_NOK) - lag(log(oil_krone$USD_NOK))

#ta vekk første raden (NA verdier)
oil_krone = oil_krone[-1,]

pacf(oil_krone$oil_price_ret)

pacf(oil_krone$USD_NOK_ret)

Nå kan vi kjøre en model. Basert på undersøkelser av ACF-figurer og tester for autokorrelasjon i restleddet, kommer jeg fram til følgende dynamisk modell:

\[\Delta k_t = \beta_1 \Delta k_{t-1} + \beta_2 \Delta k_{t-2} + \gamma_1 \Delta p_{t} + \gamma_2 \Delta p_{t-1} + \epsilon_t\]

Der \(\Delta k_t\) er avkastningen av kronekursen og \(\Delta p_t\) er avkastningen av oljeprisen.

Resultatene finner man under.

krone_mod = dynlm(USD_NOK_ret ~ lag(USD_NOK_ret) + lag(USD_NOK_ret,2) + oil_price_ret + lag(oil_price_ret), data=oil_krone)
summary(krone_mod)
## 
## Time series regression with "numeric" data:
## Start = 1, End = 213
## 
## Call:
## dynlm(formula = USD_NOK_ret ~ lag(USD_NOK_ret) + lag(USD_NOK_ret, 
##     2) + oil_price_ret + lag(oil_price_ret), data = oil_krone)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.076065 -0.014322 -0.000804  0.015159  0.063293 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         -0.0002711  0.0014860  -0.182   0.8554    
## lag(USD_NOK_ret)     0.3168103  0.0683603   4.634 6.31e-06 ***
## lag(USD_NOK_ret, 2) -0.1143092  0.0615892  -1.856   0.0649 .  
## oil_price_ret        0.1375087  0.0176926   7.772 3.50e-13 ***
## lag(oil_price_ret)  -0.0284003  0.0192089  -1.478   0.1408    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.02164 on 208 degrees of freedom
##   (2 observations deleted due to missingness)
## Multiple R-squared:  0.3362, Adjusted R-squared:  0.3235 
## F-statistic: 26.34 on 4 and 208 DF,  p-value: < 2.2e-16

Noen viktige ting å se her:

  • Valutakursen er sterkt påvirket av sin lag fra en måned siden og to måneder siden, men i ulike retninger. Dette betyr at hvis kronekursen øker en måned, så kan vi forvente at den øker neste måned også, men også at den kommer litt tilbake etter det.

  • En endring i oljeprisen har en signifikant og positiv effekt på valutakursen med dollaren. Vi kan tolke det sånt: “en 10% endring (økning) i oljeprisen fører til, i gjennomsnitt, en .8% endring i valutakursen.”

  • Adjusted R-squared er 0,23. Vi kan tolke dette med å si at vi har forklart 23% av variasjonen i valutakursen med vår enkel modell. Men vær forsiktig med å være for opptatt av ulike R-squared statistikk. Målet med regresjonsanalyser er sjeldent å få høyest R-squared. Vi kan ha en god model som gir innsikt i et økonomisk forhold, uten at modellen har særlig høy R-squared.