Labb 9 Introduksjon til tidsrekke-regresjon

Makroøkonomi handler ofte om endringer over tid. Vi har lyst å vite hvorfor BNP, arbeidsledighet, inflasjon, osv, faller eller øker over tid. Vi er intressert i hvordan faktorer som pengepolitikk, finanspolitikk og økonomiens struktur påvirker disse variablene på kort-, medium- og langsikt. Da er vi avhengig av en spesialisert form av regresjonsanalyse som tar tid og dynamikk i betrakning. For å være enda mer spesifikk - vi må ta i betrakning at våre variabler er organisert i en tidsrekke, og at dette kan spille en stor rolle i hvordan vi kjører og tolker regresjoner.

Følgende introduksjon er delvis basert på Modern Econometrics av Wooldridge.

Distribuerte lag modeller

Når vi har en y-variabel som er modellert på en x-variabel og lagget x-variabler, så kaller vi det for en distribuerte lag modeller. Det kan skrives:

\(y_t = \alpha + \delta_0 x_t + \delta_1 x_{t-1} + \delta_2 x_{t-2} + u\)

La oss si at vi vil analysere en midlertidig effekt: en “sjokk” som deretter avtar i effekt. Vi kan tenke på, for eksempel, en hopp i oljeprisen påvirker inflasjon over tid.

• Vi sier at sjokke treffer på tid \(t\). I tidsperioden før t, så er x variablen konstant, \(x=c\).

Vi begynner med perioden før sjokket treffer:

\(y_{t-1} = \alpha_0 + \delta_0 c + \delta_1 c + \delta_2 c\)

I periode t får vi:

\(y_t = \alpha_0 + \delta_0 (c+1) + \delta_1 c + \delta_2 c\)

og så videre:

\(y_{t+1} = \alpha_0 + \delta_0 c + \delta_1 (c+1) + \delta_2 c\)

\(y_{t+2} = \alpha_{0} + \delta_{0} c + \delta_{1} c + \delta_{2} (c+1)\)

\(y_{t+3} = \alpha_0 + \delta_0 c + \delta_1 c + \delta_2 c\)

Vi kan da se hva effekten av sjokken er over tid, det blir kalt “Impact Propensity”

tid	impact propensity
t	\(y_t-y_{t-1} = \delta_0\)
t+1	\(y_{t+1}-y_{t-1} = \delta_1\)
t+2	\(y_{t+2} - y_{t-1} = \delta_2\)
t+3	\(y_{t+3} - y_{t-1} = 0\)

Hva om det er en permanent effekt?

Igjen, i periode t får vi:

\(y_t = \alpha_0 + \delta_0 (c+1) + \delta_1 c + \delta_2 c\)

og så videre:

\(y_{t+1} = \alpha_0 + \delta_0 (c+1) + \delta_1 (c+1) + \delta_2 c\)

\(y_{t+2} = \alpha_0 + \delta_0 (c+1) + \delta_1 (c+1) + \delta_2 (c+1)\)

\(y_{t+3} = \alpha_0 + \delta_0 (c+1) + \delta_1 (c+1) + \delta_2 (c+1)\)

Da snakker vi om “Long run propensity”: permanent effekt.

tid	long-run propensity
t	\(y_t-y_{t-1} = \delta_0\)
t+1	\(y_{t+1}-y_{t-1} = \delta_0 + \delta_1\)
t+2	\(y_{t+2} - y_{t-1} = \delta_0 + \delta_1 + \delta_2\)
t+3	\(y_{t+3} - y_{t-1} = \delta_0 + \delta_1 + \delta_2\)

Seriekorrelasjon

Hvis vi vil ha korrekt standard feil i våre regresjoner, så kan vi ikke har seriekorrelasjon i våre feilledd.

Seriekorrelasjon defineres som:

\(corr(u_t, u_s | X)\)

der \(t\neq s\)

Det vil si, at det ikke er noe korrelasjon over tid i feilleddene. Intuitivt, hvis vi ser på residualene av en regresjon, så burde de se tilfeldig ut. Det burde ikke være noe korrelasjon.

Log transformasjon

I økonometri vil man ofte se at man log-transformerer data. * Man bruker log-transformasjon til å gjøre data mer linear * Hvis man tranformerer begge sidene av regresjonen, så kan man tolke resultatene som elastitet

Si at jeg vil kjøre en regresjon av BNP på pengeetterspørsel, M. Jeg log-transformerer min regresjon:

\(log(M_t) = \alpha_0 + \delta_0 * log(GDP_t) + \delta_1*log(GDP_{t-1}) + \delta_2*log(GDP_{t-2})\)

\(\delta_0\) tolkes som “kort-siktig elasitet”: \(\delta_0 = \frac{\% \Delta M}{\% \Delta GDP_t}\)

Lang-siktig elasitet kan tolkes som: \(\delta_0 + \delta_1 + \delta_2\)

• Husk: for å gjøre en log transformasjon, må alle verdiene være > 0.

Dummy-variabler i tidsrekke regresjon.

Vi kan inkludere dummy-variabler i vår tidsrekke regresjon. For eksempel, la oss si at vi vil prøve å estimere effekten av regjering på BNP. Kanskje vi kjører følgende regresjon:

\(log(BNP_t) = \beta X + \delta D_t + u_t\)

der \(D_t = [0,0,0,0,1,1,1,1,0,0,0]\)

Trender

Ofte vil vi modellere en tidstrend i våre regresjoner.

For en linear trend, kan vi skrive det som:

\(y_t = \alpha + \beta*t + u_t\)

u_t = rnorm(100,0,1) #100 tilfeldig trekk fra standard normal fordeling
alpha = 2.5
beta = .3
t=seq(1:100)
y = alpha + beta*t+u_t

plot(t,y, type="line")

## Warning in plot.xy(xy, type, ...): plot type 'line' will be truncated to
## first character

Hvis vi modellerer log(y_t), så blir den en eksponensiel trend:

\(log(y_t) = \alpha + \beta t 0 + u_t\)

alpha = 2.5
beta = .03
y= exp(alpha + beta*t) + u_t
plot(t, y, type="line")

## Warning in plot.xy(xy, type, ...): plot type 'line' will be truncated to
## first character

Tidsrekke med R

Vi laster inn noen pakker som vi trenger og en tidsrekke datasett av renter, inflasjon og statlig budsjettunderskudd:

library(foreign) #for å importere data i stata format (.dta)
library(ggplot2)
library(reshape2) #for melt
library(zoo)

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

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

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

##   year   i3  inf  rec  out        def i3_1 inf_1      def_1        ci3
## 1 1948 1.04  8.1 16.2 11.6 -4.6000004   NA    NA         NA         NA
## 2 1949 1.10 -1.2 14.5 14.3 -0.1999998 1.04   8.1 -4.6000004 0.06000006
## 3 1950 1.22  1.3 14.4 15.6  1.2000008 1.10  -1.2 -0.1999998 0.12000000
## 4 1951 1.55  7.9 16.1 14.2 -1.9000006 1.22   1.3  1.2000008 0.32999992
## 5 1952 1.77  1.9 19.0 19.4  0.3999996 1.55   7.9 -1.9000006 0.22000003
## 6 1953 1.93  0.8 18.7 20.4  1.6999989 1.77   1.9  0.3999996 0.15999997
##   cinf      cdef y77
## 1   NA        NA   0
## 2 -9.3  4.400001   0
## 3  2.5  1.400001   0
## 4  6.6 -3.100001   0
## 5 -6.0  2.300000   0
## 6 -1.1  1.299999   0

Her er i3 renter, inf inflasjon og def budsjett-underskudd. Og vi kan se at de er alle organisert med år.

Som vanlig kan det være greit å plotte datasettet:

intdef_l <- intdef[c("year","i3","inf", "def")]
intdef_l <- melt(intdef_l, id.var="year")

ggplot(intdef_l, aes(year, value)) +
  geom_line() +
  facet_wrap(~variable, nrow=3)

ggplot(intdef_l, aes(x=year, y=value, linetype=variable)) +
  geom_line() + theme_bw() + ylab("%")

Kan man si noe om forholdet mellom disse seriene ut fra plottet?

Svar

Det ser ut som at variablene følger hverandre til en viss grad. Vi vet at det burde være et forhold mellom nominelle renten og inflasjon. Vi kan kanskje også lese at det er viss forhold mellom lavere budsjett-underskudd og lavere renter.

Vi kan kjøre en vanlig regresjon og se hva resultatene blir:

model1 <- lm(i3 ~ inf + def, data=intdef)
summary(model1)

## 
## Call:
## lm(formula = i3 ~ inf + def, data = intdef)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.9948 -1.1694  0.1959  0.9602  4.7224 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.73327    0.43197   4.012  0.00019 ***
## inf          0.60587    0.08213   7.376 1.12e-09 ***
## def          0.51306    0.11838   4.334 6.57e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.843 on 53 degrees of freedom
## Multiple R-squared:  0.6021, Adjusted R-squared:  0.5871 
## F-statistic: 40.09 on 2 and 53 DF,  p-value: 2.483e-11

Hvordan kan du tolke resultatene? Finnes det problemer med denne regresjonen / tolkningen?

Svar

Både inflasjon og budsjett-underskudd er signifikant. Vi kunne tolke koefisienten på inflasjon sånn at en 1-prosentpoeng økning i inflasjon, fører til renter som er 0,6 prosentpoeng høyere, gitt budsjettunderskuddet. Vi kunne gjøre en tolkning av budsjett-underskuddet på samme vis.

Regresjonen har noen mulige problemer, som vi skal snakke om mer senere. 1.) Vi må ta hensyn til autokorrelasjon - det vil si at inflasjon i periode t er avhengig av inflasjon i peridoe t-1 2.) Kanskje viktigere her er at vi burde ta hensyn til en felles trend i datasettet.

Nå vil vi jobbe med å sette inn datasettet i riktig formatt for tidsrekke-analyse.

I denne eksempelen har vi data for hvert år. Når vi har data en gang per periode så kaller vi det “regular time series”. Da kan vi bruke r’s innebygget ts formatt

i3_ts <- ts(intdef$i3, start=1948, frequency=1)
str(i3_ts)

##  Time-Series [1:56] from 1948 to 2003: 1.04 1.1 1.22 1.55 1.77 ...

Her skrev vi at vår tidsrekke skulle begynne i 1948 og at vi har årlig data (frequency=1). Om det var månedsdata så kunne vi for eksempel skrive frequency=12 og vi kunne si hvilken måned dataen skal started ved å skrive start=c(1948,2), som betyr at det første datapunktet var februar, 1948.

ts-formattet har også innebygd plotte-funksjoner:

plot(i3_ts)

Men hva om dataene ikke er regulære. Tenk, for eksempel, på daglig aksjedata der man ikke har data i helgene eller der man generelt sett ikke har data i regelmessige intervaller. Da kan vi bruke pakken og formattet zoo

library(zoo)
intdef_zoo <- zoo(intdef, order.by = intdef$year)
plot(intdef_zoo$i3)

Når vi jobber med tidsrekke data og regresjon vil vi ofte endre og transformere dataen - lags, differensiering, tid-trender, og sesong-justering. Vi bruker pakken dynlm for dette:

install.packages("dynlm")

library(dynlm)

Distribuerte lag modeller:

Nå skal vi begynne å kjøre noen tidsrekke-regresjoner. I tillegg til dylm og zoo, vil vi bruke et par andre pakker:

library(lmtest) #for å kjøre lm-test (en type F-test)
library(car)

## Loading required package: carData

Vi bruker en datasett om fruktbarhet

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

##     gfr    pe year t tsq  pe_1 pe_2 pe_3 pe_4 pill ww2 tcu      cgfr   cpe
## 1 124.7  0.00 1913 1   1    NA   NA   NA   NA    0   0   1        NA    NA
## 2 126.6  0.00 1914 2   4  0.00   NA   NA   NA    0   0   8  1.900002  0.00
## 3 125.0  0.00 1915 3   9  0.00    0   NA   NA    0   0  27 -1.599998  0.00
## 4 123.4  0.00 1916 4  16  0.00    0    0   NA    0   0  64 -1.599998  0.00
## 5 121.0 19.27 1917 5  25  0.00    0    0    0    0   0 125 -2.400002 19.27
## 6 119.8 23.94 1918 6  36 19.27    0    0    0    0   0 216 -1.199997  4.67
##   cpe_1 cpe_2 cpe_3 cpe_4 gfr_1    cgfr_1    cgfr_2    cgfr_3   cgfr_4
## 1    NA    NA    NA    NA    NA        NA        NA        NA       NA
## 2    NA    NA    NA    NA 124.7        NA        NA        NA       NA
## 3  0.00    NA    NA    NA 126.6  1.900002        NA        NA       NA
## 4  0.00     0    NA    NA 125.0 -1.599998  1.900002        NA       NA
## 5  0.00     0     0    NA 123.4 -1.599998 -1.599998  1.900002       NA
## 6 19.27     0     0     0 121.0 -2.400002 -1.599998 -1.599998 1.900002
##   gfr_2
## 1    NA
## 2    NA
## 3 124.7
## 4 126.6
## 5 125.0
## 6 123.4

Datasettet ser på gjennomsnittlig fruktbarhet (gfr) over tid i USA. Det aktuelle spørsmålet er om skattefradrag for barn (pe) vil ha en effekt. Man kan kontrollere for andre verdenskrig (ww2) og for oppfinnelsen av P-piller (pill).

Vi setter datasettet i zoo-formatt:

zoo_fertil3 <- zoo(fertil3, order.by = fertil3$year)
plot(zoo_fertil3$gfr)

Nå kan vi prøve å kjøre en dynamisk modell.

ts_model1 <- dynlm(gfr ~ pe + L(pe) + L(pe,2) + ww2 + pill, data=zoo_fertil3) 
summary(ts_model1)

## 
## Time series regression with "zoo" data:
## Start = 1915, End = 1984
## 
## Call:
## dynlm(formula = gfr ~ pe + L(pe) + L(pe, 2) + ww2 + pill, data = zoo_fertil3)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -24.6461  -9.5409  -0.0312   8.3378  29.1295 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  95.87050    3.28196  29.211  < 2e-16 ***
## pe            0.07267    0.12553   0.579   0.5647    
## L(pe)        -0.00578    0.15566  -0.037   0.9705    
## L(pe, 2)      0.03383    0.12626   0.268   0.7896    
## ww2         -22.12650   10.73197  -2.062   0.0433 *  
## pill        -31.30499    3.98156  -7.862 5.63e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 14.27 on 64 degrees of freedom
## Multiple R-squared:  0.4986, Adjusted R-squared:  0.4594 
## F-statistic: 12.73 on 5 and 64 DF,  p-value: 1.353e-08

Dette er som en standard linear modell, som vi estimerer med OLS og funksjonen lm, men vi tar også i betrakning at effekten av skattefradrag kan være litt forsinket, så derfor inkluderer vi leddene L(pe) og L(pe, 2) som betyr at man skal bruke effekten av skattefradraget fra året før og 2 år før (\(pe_{t-1}, pe_{t-2}\)) i regresjonen.

Kan du tolke resultatene? Ser det ut som skattefradrag har noen effekt?

P-verdiene sier kun noe om de individuelle leddene er signifikante. La oss teste om alle leddene med pe, sammen, kan vises å ha en effekt.

#modell med restriksjoner (pe, pe_t-1, pe_t-2 == 0 )
ts_model1_res <- dynlm(gfr ~ ww2 + pill, data=zoo_fertil3) 

#F-test
r2_res <- summary(ts_model1_res)$r.squared
r2_ures <- summary(ts_model1)$r.squared
n <- length(zoo_fertil3$gfr)
k <- 5 #antall variabler på høyre siden, unrestricted
r <- 3 #antall restriksjoner
f_stat <- (r2_ures - r2_res)/(1-r2_ures)*(n-k-1)/r
f_stat

## [1] 3.738784

Da kan vi finne en p-verdi for F-stat ved bruk av F-fordelingen

Vi kan finne p-verdien av en f-statistikk med r-restriksjoner og n-k-1 frihetsgrad:

1-pf(f_stat, r, (n-k-1))

## [1] 0.01520076

Vi kan visualisere hva det betyr med en F-fordeling of F-test. Vi begynner ved å simulere 100.000 trekk fra en F-fordeling med r restriksjoner og (n-k-1) frihetsgrad:

f_sim<-data.frame(f_sim = rf(100000, r, (n-k-1)))
head(f_sim)

##       f_sim
## 1 0.2941641
## 2 0.3133305
## 3 0.4274549
## 4 0.9555071
## 5 1.2844487
## 6 4.5248809

Vi kan da simulere fordelingen ved å vise en histogram av alle de simulerte trekkene fra fordelingen.

ggplot(f_sim) + 
  geom_histogram(aes(f_sim), bins=50) +
  geom_vline(xintercept=f_stat, color="red")

length(f_sim$f_sim[f_sim>f_stat])/length(f_sim$f_sim)*100

## [1] 1.496

Vi kan også bruke kommandoen linearHypothesis fra car-pakken til å utføre en F-test med bare en linje:

linearHypothesis(ts_model1, c("pe + L(pe) + L(pe, 2) = 0")) 

## Linear hypothesis test
## 
## Hypothesis:
## pe  + L(pe)  + L(pe, 2) = 0
## 
## Model 1: restricted model
## Model 2: gfr ~ pe + L(pe) + L(pe, 2) + ww2 + pill
## 
##   Res.Df   RSS Df Sum of Sq      F   Pr(>F)   
## 1     65 15358                                
## 2     64 13033  1    2325.8 11.421 0.001241 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Vi skriver inn tsmodel1 som modellen uten restriksjoner, og pe + L(pe) + L(pe, 2)=0 som modellen med restriksjoner.

Du kan lese resultatene som at modellen med pe og laggene av pe gjør en bedre jobb av å forklare fruktbarhet en den enklere modellen.

Trends

Huspriser har en tendens til å øke over tid, mens renten og inflasjon har hatt en tendens å synke gjennom de siste 30 årene. Disse trendene burde vi kontrollere for når vi utfører tidsrekke-regresjon.

Vi begynner med å laste inn data om investering i boliger og boligpriser fra USA. Vi skal også introdusere en pakke som heter stargazer som lar oss automatisk lage fine tabeller fra våre resultater.

install.packages("stargazer")

Vi laster in vår data om boliginvestering i USA.

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

##   year   inv    pop  price     linv     lpop     lprice t     invpc
## 1 1947 54864 144126 0.8190 10.91261 11.87844 -0.1996712 1 0.3806669
## 2 1948 64717 146631 0.8649 11.07778 11.89567 -0.1451414 2 0.4413596
## 3 1949 63150 149188 0.8456 11.05327 11.91296 -0.1677088 3 0.4232914
## 4 1950 86014 151684 0.8765 11.36227 11.92955 -0.1318186 4 0.5670605
## 5 1951 70610 154287 0.8819 11.16493 11.94657 -0.1256766 5 0.4576536
## 6 1952 68574 156954 0.8842 11.13567 11.96371 -0.1230720 6 0.4369051
##       linvpc   lprice_1   linvpc_1     gprice     ginvpc
## 1 -0.9658305         NA         NA         NA         NA
## 2 -0.8178953 -0.1996712 -0.9658305  0.0545298  0.1479352
## 3 -0.8596944 -0.1451414 -0.8178953 -0.0225674 -0.0417991
## 4 -0.5672894 -0.1677088 -0.8596944  0.0358903  0.2924051
## 5 -0.7816427 -0.1318186 -0.5672894  0.0061420 -0.2143534
## 6 -0.8280393 -0.1256766 -0.7816427  0.0026046 -0.0463966

Nå vil vi begrense vår datasett til kun variablene “year”, “invpc” (boliginvestering per kapita), “price” (prisindek for boliger). Vi bruker funksjonen melt til å få datasetten i lang-formatt der alle variabel-navnene ligger i en kolonne. Dette er formattet vi trenger for å bruke ggplot2.

hseinv_l <- hseinv[c("year", "invpc", "price")]
hseinv_l <- melt(hseinv_l, id.var="year")

ggplot(hseinv_l, aes(year, value)) +
  geom_line() +
  facet_wrap(~variable, nrow=2, scales="free_y")

Vi vil nå omformatterer det opprinnelige datasettet (i bred-formatt) til tidsrekke med zoo:

hseinv_zoo <- zoo(hseinv, order.by=hseinv$year)

Nå kan vi kjører noen dynamiske modeller:

bolig_model1 <- dynlm(log(invpc) ~ log(price), data =hseinv_zoo)
bolig_model2 <- dynlm(log(invpc)~log(price) + trend(hseinv_zoo), data=hseinv_zoo)
summary(bolig_model1)

## 
## Time series regression with "zoo" data:
## Start = 1947, End = 1988
## 
## Call:
## dynlm(formula = log(invpc) ~ log(price), data = hseinv_zoo)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.45991 -0.08694 -0.01264  0.08651  0.34672 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.55023    0.04303 -12.788 1.03e-15 ***
## log(price)   1.24094    0.38242   3.245  0.00238 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1554 on 40 degrees of freedom
## Multiple R-squared:  0.2084, Adjusted R-squared:  0.1886 
## F-statistic: 10.53 on 1 and 40 DF,  p-value: 0.002376

summary(bolig_model2)

## 
## Time series regression with "zoo" data:
## Start = 1947, End = 1988
## 
## Call:
## dynlm(formula = log(invpc) ~ log(price) + trend(hseinv_zoo), 
##     data = hseinv_zoo)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.45092 -0.08583 -0.01734  0.08517  0.26024 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       -0.913060   0.135613  -6.733    5e-08 ***
## log(price)        -0.380961   0.678835  -0.561  0.57787    
## trend(hseinv_zoo)  0.009829   0.003512   2.798  0.00794 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1436 on 39 degrees of freedom
## Multiple R-squared:  0.3408, Adjusted R-squared:  0.307 
## F-statistic: 10.08 on 2 and 39 DF,  p-value: 0.000296

Hvordan kan man tolke disse resultatene?

Svar

I første regresjonen er pris signifikant. Men når vi inkluderer en trend, så er det ikke lenger signifkant.

Hvorfor ville det være feil å ekskludere en trend i dette tilfellet?

Svar

Det kan være en uobservert variabel som påvirker begge variablene over tid. Dette vil vi kontrollere for.

Sesongjustering

Vi vet at boligpriser er altid litt dyrere i august og at energietterspørsel er høyere om vinteren (i norge) når alle skal varme opp husene. Disse sesong-effekter burde man ofte prøve å kontrollere.

Generelt sett, når vi snakker om sesongjustering, snakker vi om alle effekter som går i en viss syklus - årlig, månedlig, daglig, osv. Det finnes faktisk en hel bransje av statistikk som handler om hvordan å håndere sesongeffekter.

Som eksempel, vi bruker data på internasjonal-handel av metallen barium.

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

##     chnimp   bchlimp befile6 affile6 afdec6 befile12 affile12 afdec12
## 1 220.4620  9578.376       0       0      0        0        0       0
## 2  94.7980 11219.480       0       0      0        0        0       0
## 3 219.3575  9719.900       0       0      0        0        0       0
## 4 317.4215 12920.950       0       0      0        0        0       0
## 5 114.6390  9790.446       0       0      0        0        0       0
## 6 129.5240 11020.470       0       0      0        0        0       0
##   chempi        gas rtwex spr sum fall  lchnimp     lgas   lrtwex  lchempi
## 1  100.1 7830000128 86.74   0   0    0 5.395725 22.78123 4.462915 4.606170
## 2  100.9 8819999744 85.63   1   0    0 4.551748 22.90029 4.450036 4.614130
## 3  101.1 8449999872 85.42   1   0    0 5.390703 22.85743 4.447580 4.616110
## 4  102.5 9240000512 87.29   1   0    0 5.760231 22.94681 4.469236 4.629863
## 5  104.1 9150000128 86.60   0   1    0 4.741788 22.93702 4.461300 4.645352
## 6  104.8 9520000000 84.63   0   1    0 4.863866 22.97666 4.438289 4.652054
##   t feb mar apr may jun jul aug sep oct nov dec   percchn
## 1 1   1   0   0   0   0   0   0   0   0   0   0 2.3016641
## 2 2   0   1   0   0   0   0   0   0   0   0   0 0.8449411
## 3 3   0   0   1   0   0   0   0   0   0   0   0 2.2567880
## 4 4   0   0   0   1   0   0   0   0   0   0   0 2.4566419
## 5 5   0   0   0   0   1   0   0   0   0   0   0 1.1709270
## 6 6   0   0   0   0   0   1   0   0   0   0   0 1.1753041

Her kan man se at datasettet inkluderer måneds-dummyvariabler. Da er det enkelt å kontrollere sesongeffekten

Vi begynner med å sette det i tidsrekke-formatt.

barium_ts <- ts(barium, start= c(1978,2), frequency=12)

barium_model <- dynlm(log(chnimp) ~ log(chempi) + log(gas) + log(rtwex) + befile6 + affile6 + afdec6 + season(barium_ts), data=barium_ts) 
summary(barium_model)

## 
## Time series regression with "ts" data:
## Start = 1978(2), End = 1988(12)
## 
## Call:
## dynlm(formula = log(chnimp) ~ log(chempi) + log(gas) + log(rtwex) + 
##     befile6 + affile6 + afdec6 + season(barium_ts), data = barium_ts)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.98535 -0.36207  0.07366  0.41786  1.37734 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          16.779215  32.428645   0.517   0.6059    
## log(chempi)           3.265062   0.492930   6.624 1.24e-09 ***
## log(gas)             -1.278140   1.389008  -0.920   0.3594    
## log(rtwex)            0.663045   0.471304   1.407   0.1622    
## befile6               0.139703   0.266808   0.524   0.6016    
## affile6               0.012632   0.278687   0.045   0.9639    
## afdec6               -0.521300   0.301950  -1.726   0.0870 .  
## season(barium_ts)Feb -0.417711   0.304444  -1.372   0.1728    
## season(barium_ts)Mar  0.059052   0.264731   0.223   0.8239    
## season(barium_ts)Apr -0.451483   0.268386  -1.682   0.0953 .  
## season(barium_ts)May  0.033309   0.269242   0.124   0.9018    
## season(barium_ts)Jun -0.206332   0.269252  -0.766   0.4451    
## season(barium_ts)Jul  0.003837   0.278767   0.014   0.9890    
## season(barium_ts)Aug -0.157064   0.277993  -0.565   0.5732    
## season(barium_ts)Sep -0.134161   0.267656  -0.501   0.6172    
## season(barium_ts)Oct  0.051693   0.266851   0.194   0.8467    
## season(barium_ts)Nov -0.246260   0.262827  -0.937   0.3508    
## season(barium_ts)Dec  0.132838   0.271423   0.489   0.6255    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.6012 on 113 degrees of freedom
## Multiple R-squared:  0.3583, Adjusted R-squared:  0.2618 
## F-statistic: 3.712 on 17 and 113 DF,  p-value: 1.282e-05

Du kan notere deg at dette er akkurat likt å bare inkludere dummy-variablene i selve regresjonen.

Oppgave: Pengetilbud og priser

Vi skal laste ned data fra Norges Bank på historisk pengemengde i økonomien og inflasjon (KPI).

Her laster vi inn datasettene fra kursiden.

penger = read.csv("http://jmaurit.github.io/anvendt_macro/data/money.csv", sep=";", stringsAsFactors = FALSE)
priser = read.csv("http://jmaurit.github.io/anvendt_macro/data/KPI.csv", sep=";", dec=",", stringsAsFactors = FALSE)

Vi formatterer til date, og drar ut KPI serien. Vi ser kun på data etter 1950.

priser["date"] = as.Date(priser$Year, format="%d.%m.%Y")
KPI = priser[c("date", "CPI")]
KPI = KPI[KPI$date>=as.Date("1950-01-01"),]
head(KPI)

##            date  CPI
## 2077 1950-01-31 5.05
## 2078 1950-02-28 5.05
## 2079 1950-03-31 5.08
## 2080 1950-04-30 5.22
## 2081 1950-05-31 5.26
## 2082 1950-06-30 5.26

Vi formatterer også en tilsvarende date variable for pengemengde, og drar ut M2_broad. Dette er en “bred” definisjon av penger: det inkluderer kontanter og brukskonto.

penger["date"] = as.Date(penger$Year, format="%d.%m.%Y")

M2 = penger[c("date", "M2_broad")]
M2 = M2[M2$date>=as.Date("1950-01-01"),]

Her skal vi bruke merge til å sette sammen dataseriene. Vi formaterer date kolonen til å være i standard år-måned formatt. Vi lagrer dataframen lokalt som en .csv fil med write.csv

likv = merge(M2, KPI, by="date")
likv$date   = as.Date(as.yearmon(likv$date))

write.csv(likv, file="likv.csv", row.names=FALSE)

• Plotte de to seriene sammen

• Skap nye dataframe i ts-formatt som heter likv_ts

• Transformere seriene først til log og da endring (per år). Plotte de to seriene først som tidsrekker og da mot hverandre. Hva ser man? Hint:
• Bruk diff til å skape serie som er differanser (endring). Bruke lag=12, hvis man vil sammenligne endring fra samme måneden året før.

• For å bruke ggplot, må man først omformattere ts-dataframen til data.frame-formatt)

• Kjør flere dynamisk regresjon av M2 på KPI med ulike lagger av KPI. Tolk resultatene. Prøv å inkludere en trend i regresjonen.

• Kjør regresjoner der venstre siden er endring i log CPI over et år (bruk diff()). Tolk resultatene.

• Kjør regresjoner der også høyre-side variablen (M2_broad) er i årlig endring-form. Tolk resultatene.