Persistens og stasjonæritet

Intro:

Når vi bruker tidsrekke-data, så bryter vi en av hoved-antakelsene som vi vanligvis bruker til å kjøre en regresjonsanalyse: tilfeldig utvalg.

Vi husker at tilfeldig utvalg sier noe om at dataens rekkefølgen ikke skal spille noe rolle. Vi kunne tatt vår data, kjørt en regresjon, deretter kunne vi byttet rekkefølgen, og kjørt samme regresjon og fått omtrent samme resultat.

Med tidsrekke-data, så er tilfeldig utvalg åpenbart noe som ikke holder. Rekkefølgen–tidsperioden–står helt sentralt i analysen. Heldigvis kan vi fortsatt bruke regresjon med tidsrekke data. Men vi må være forsiktig med å sjekke visse betingelser og eventuelt gjøre noen transformasjoner.

Stasjonæritet er den viktigiste av disse betingelsene. Intuitivt, sier stasjonæritet at vår tidsrekke ser ut som at det ligner på tilfeldig utvalg: at vi kunne byttet rekkefølgen av vår datasettet, og den ville fortsatt sett ganske likt ut. For eksempel, stasjonær data har ikke en tidstrend.

Den viktigiste tingen du burde vite om stasjonæritet (som vi kommer til å definere og forklare) er at hvis din data er ikke-stasjonær, så kan du ikke stole på resultatene fra en regresjon uten videre.

Det er verdt å si det en gang til:

Hvis din data er ikke-stasjonær, så kan du ikke stole på resultatene fra en regresjon.

Intro til stasjonæritet og persistens

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.numeric
library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.0 ──
## ✓ ggplot2 3.3.0     ✓ purrr   0.3.4
## ✓ tibble  3.0.3     ✓ dplyr   0.8.5
## ✓ tidyr   1.1.0     ✓ 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(zoo)
library(fpp2)
## Loading required package: forecast
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
## Loading required package: fma
## Loading required package: expsmooth
library(stargazer)
## 
## 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=stargazer

Vi starter med å se på et begrep som er relatert til stasjonæritet: persistens.

Vi bruker data om gjennomsnitts-priser og -avkastninger på New York Stock Exchange. Dette er data i Stata-formatt som man kan laste ned med linken: http://fmwww.bc.edu/ec-p/data/wooldridge/nyse.dta.

Her er hva de første radene av vår datasett ser ut som. Her ser vi at vi har både daglig priser og avkastninger (“returns” - % endringer fra dag til dag). Datasettet har også lagget verdier - det vil si at hele kolonen er skiftet nedover.

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

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.9600030

Nå kan vi kjøre en enkel AR regresjon. Det man tester med denne regresjonen er om avkastning på New York Stock Exchange en dag kan si noe om avkastningen neste dag og dagen etter det, osv. Dette er definisjonen av persistens. Hvis en serie ikke er persistent, så er det lite eller ingen sammenheng - avkastningen en dag er uavhengig av avkastningen dagen før og dagen før det, osv.

\[ R_t = \alpha + \beta_1 R_{t-1} + \beta_2 R_{t-2} + ... + \epsilon\]

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.01

Resultatene fra en rekke modeller med ulike antall lagger viser at det er en forholdsvis sterk forhold mellom avkastningen en dag og den neste. Men etter det så blir det fort lite sammenheng. Det er, med andre ord, persistens som varer omtrent en dag i avkastningene.

Persistens i seg selv er ikke et stort problem (se materialene for hvordan å håndtere autokorrelasjon.) Så lenge den persistensen dør ut ganske fort, fra en periode til den neste.

Problemet er hvis det blir veldig mye persistens (“highly persistent”), så at en endring i en periode har en effekt på størrelsene en lang tid i etterkant.

La oss se på en annen tidsrekke: her er pris-data for GE-aksjen i USA

ge_stock = read_csv("https://jmaurit.github.io/anvendt_macro/data/GE.csv")
## Parsed with column specification:
## cols(
##   Date = col_date(format = ""),
##   Open = col_double(),
##   High = col_double(),
##   Low = col_double(),
##   Close = col_double(),
##   `Adj Close` = col_double(),
##   Volume = col_double()
## )
head(ge_stock)
## # A tibble: 6 x 7
##   Date        Open  High   Low Close `Adj Close`   Volume
##   <date>     <dbl> <dbl> <dbl> <dbl>       <dbl>    <dbl>
## 1 2012-09-17  21.9  22.0  21.9  22.0        18.5 79283200
## 2 2012-09-18  22.0  22.2  22.0  22.2        18.7 38831200
## 3 2012-09-19  22.3  22.5  22.2  22.4        18.8 40462600
## 4 2012-09-20  22.2  22.5  22.1  22.4        19.0 43677000
## 5 2012-09-21  22.5  22.7  22.5  22.5        19.1 66551700
## 6 2012-09-24  22.4  22.5  22.3  22.4        18.9 36799800
ge_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. Tidsrekken er preget av høy auto-korrelasjon og variansen endrer seg over tid. En slik tidsrekke er ikke-stasjonær.

Hvis vi prøvde å bytte på rekkefølgen av denne tidsrekken, så ville figuren sett helt annerledes ut.

Dette gjør det sånt at man ikke kan bruke vanlig regresjon (MLM) på en konsistent måte uten videre.

Men hvis vi transformerer datarekken til daglig avkastninger - det vil si prosent endring per dag, så får vi en datarekke som ser sånt ut:

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)

Med en transformasjon av dataen så har vi nå en tidsrekke som ser mye mindre peristent ut. Her har vi nå en stasjonær tidsrekke som vi kan modellere med regresjon.

Random Walk

Den enkleste formen av en veldig persistent (“highly persistent”) tidsrekke er en random walk. Du kan lese mer her. 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 ofte brukt til å modellere oljepriser, aksjepriser, og også noen makroøkonomiske variabler.

#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 = as_tibble(y)
ycols = colnames(y)
y["t"] = 1:50
y_l = gather(y, ycols, key="variable", value="value")
## Note: Using an external vector in selections is ambiguous.
## ℹ Use `all_of(ycols)` instead of `ycols` to silence this message.
## ℹ See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.
## This message is displayed once per session.
ggplot(y_l, aes(x=t, y=value, group=variable)) +
  geom_line(alpha=.5)

I figuren vises 30 simulerte tidsrekker av 50 perioder. Man kan se at variansen øker over tid. De individuelle tidsrekkene beveger seg ikke mot en bestemt gjennomsnittsverdi, men kan heller ende opp både høyt oppe og langt nede. Seriene er veldig persistent der rekkefølgen er avgjørende. Her har vi ikke-stasjonære tidsrekker.

En annen form av random walk er random walk med drift. Det vil si at random walk mønsteret har en bestemt retning - en tidstrend. 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 = as_tibble(y)
ycols = colnames(y)
y["t"] = 1:50
y_l=gather(y, ycols, key="variable", value="value")

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 lang 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 /ikke-stasjonæritet

Highly persistent data/ikke-stasjonær data (random walk, osv) er, som nevnt, vanskelig å jobbe med i regresjons-sammenheng. Heldigvis er det ofte en enkel løsning som vi allerede har sett: Differensiering (“first difference)”.

Du kan lese mer her.

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:

y=matrix(,nrow=50, ncol=30)
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 = as.tibble(y)
ycols = colnames(y)
y["t"]=1:50
yz = zoo(y, order.by = y$t)
Dy = diff(yz)

Dy = as.tibble(Dy)
Dy["t"] = 1:49

Dy_l=Dy %>% gather(ycols, key="variable", value="value")

ggplot(Dy_l, aes(x=t, y=value, group=variable)) +
  geom_line(alpha=.2) +
  geom_hline(yintercept=2)

Nå kan vi se at variansen til tidsrekken ser omtrent konstant gjennom tid. Vi kunne klart å beregne en OLS regresjon som ville gitt oss omtrent en vannrett linje.

Stasjonæritet: Litt mer formelt

Stasjonæritet er en av de viktigste begrepene i tidsrekke-analyse. Rent intuitivt, er en stasjonær tidsrekke en tidsrekke som er stabil over tid og at rekkefølgen av data ikke spiller en avgjørende rolle. Litt mer formelt, definerer vi kovarians-stasjonæritet som:

Vi har allerede nevnt noen eksempler på tidsrekker som ikke er stasjonære:

Hvordan bruker vi stasjonæritet?

Skal vi prøve å lære noe om forholdet mellom to variabler som er meningsfylt, burde forholde mellom disse variablene være stabile over tid.

Kriterier for stasjonæritet

I vanlig regresjon krevde vi tilfeldig utvalg som kriterie for at regresjonen skal være meningsfylt. Vi vet at en tidsrekke vanligvis ikke er et tilfeldig utvalg (data er i en tidsrekke, dvs ikke tilfeldig rekkefølge.)

Istedenfor vil vi se etter det vi kaller “Weak dependence”, eller svak avhengighet:

\(corr(x_t, x_{t+h})\rightarrow0\) når h blir stor.

En annen måte å tenke på dette er at en sjokk til en serie vil dø ut og ikke bli permanent.

AR(1) prosess

La oss gå tilbake til en AR(1) prosess:

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

Først, setter vi \(\rho=.5\) og simulerer:

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

Nå la oss sette \(\rho=1\) og simulere:

y=(1)
rho = 1 #AR ledd
for(t in 2:100){
  y[t] = rho*y[t-1] + u[t] 
}
t=seq(1:100)
plot(t,y, type="line")

Og til slutt, setter vi \(\rho=1,05\) og simulerer:

u=rnorm(100,0,1)
y=(1)
rho = 1.05 #AR ledd
for(t in 2:100){
  y[t] = rho*y[t-1] + u[t] 
}
t=seq(1:100)
plot(t,y, type="line")

Hvilken av disse tre seriene ser stabil ut?

Generelt, kan vi si at hvis \(|\rho|<1\) så er serien stasjonær og weakly dependent (og ikke “highly persistent”).

Transformasjon av persistent data.

Hvis vi har en tidsrekke som med høy persistens (“highly persistent”), kaller vi det ofte \(I(1)\): integrated of order one, eller en unit-root process.

Det vi ønsker er \(I(0)\): Integrated of order 0, eller weakly persistent.

Vi har allerede sett måten vi kan gjøre dette: Transformere til endringer (eller hvis det er pris-data, avkastninger).

Her for eksempel simulerer vi igjen med en AR(1) modell der \(\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)
plot(t,y, type="line")

Nå tar vi vår serie: \([y_1, y_2, y_3....)\) og skaper en ny serie med endringer (“first difference”):

\[ \Delta y_1 = y_2 - y_1\\ \Delta y_2 = y_3 - y_2\\ \Delta y_3 = y_4 - y_3 \\ .\\ .\\ .\\ \] Vi plotter den nye serien med endringer:

diff_y = diff(y)
plot(t[2:100], diff_y, type="line")

Vi ser nå at serien ser mer stabil ut - vi har transformert vår ikke-stasjonær serie til en form som er stasjonær.

Men en sånn transformasjon er ikke magisk. Man må nå tolke resultater som endringer (“En endring i y fører til en endring i x”). Vi har også kvittet oss med noe av informasjonen som var i vår opprinnelig data-sett. Men når vi ønsker å bruke en ikke-stasjonær data-sett i en regresjon, er det ofte verdt kostnaden å gjøre denne transformasjonen.