Lab 15: Øving til eksamen
Johannes Mauritzen
Høsten 2019
I denne labben skal vi øve på å behandle, analysere og tolke tidsrekke-data. Innholdet i denne labben er i stor grad eksamensrelevant.
library(tidyverse)
library(foreign)
library(dynlm)
library(zoo)Oppgave 1: gjeld
Gjeld i en økonomi er en viktig makroøkonomisk indikatør. Her laster vi ned data om gjeld hos husholdninger.
#laster ned fra min webside
gjeld_HH = read_csv("http://jmaurit.github.io/norwayeconomy/data_series/household_debt.csv")## Warning: Missing column names filled in: 'X1' [1]## Parsed with column specification:
## cols(
## X1 = col_double(),
## currency = col_character(),
## sector = col_character(),
## credit_source = col_character(),
## time = col_date(format = ""),
## variable = col_character(),
## value = col_double()
## )#begrenser datasett til tre kolonner
gjeld_HH = gjeld_HH[c("time", "value", "credit_source")]
#omformatterer til Date formatt
gjeld_HH["time"] = as.Date(gjeld_HH$time)
#omdefinerer kategoriene for gjeldskilde
levels(gjeld_HH$credit_source) = c("Lån fra statlige låneinstitutter", "Lån fra banker", "Lån fra kredittforetak")
head(gjeld_HH)## # A tibble: 6 x 3
## time value credit_source
## <date> <dbl> <chr>
## 1 1987-12-01 246. L202 Lσn fra banker
## 2 1988-01-01 246. L202 Lσn fra banker
## 3 1988-02-01 249. L202 Lσn fra banker
## 4 1988-03-01 250. L202 Lσn fra banker
## 5 1988-04-01 251. L202 Lσn fra banker
## 6 1988-05-01 252. L202 Lσn fra bankerHer ser vi at husholdningers gjeld kommer hovedsakelig fra banker og kredittforetak, som inkluderer både boligkredittforetak og andre kredittforetak som tilbyr, for eksempel, forbrukslån. Mer info [her][https://www.finanstilsynet.no/forbrukerinformasjon/bank-og-finans/]
ggplot(gjeld_HH, aes(x=time, y=value, color=credit_source)) +
geom_line()
Vi vil lage en serie som summer over alle gjeldskildene. Vi bruker spread funksjonen fra pakken tidyverse
gjeld_HH_wide = pivot_wider(gjeld_HH, names_from ="credit_source", values_from = "value")
gjeld_HH_wide["total_hh_gjeld"] = rowSums(gjeld_HH_wide[, c(2:4)])
head(gjeld_HH_wide)## # A tibble: 6 x 5
## time `L202 Lσn fra ba… `L201 Lσn fra st… `L203 Lσn fra k… total_hh_gjeld
## <date> <dbl> <dbl> <dbl> <dbl>
## 1 1987-12-01 246. 97.6 40.5 384.
## 2 1988-01-01 246. 99.0 41.3 386.
## 3 1988-02-01 249. 99.7 41.6 390.
## 4 1988-03-01 250. 100. 42.3 393.
## 5 1988-04-01 251. 100. 43.8 395.
## 6 1988-05-01 252. 101. 45.1 399.Nå kan vi formattere serien som en zoo-tidsrekke:
gjeld_HH_zoo = zoo(gjeld_HH_wide$total_hh_gjeld, order.by=gjeld_HH_wide$time)
plot(gjeld_HH_zoo)
Hva sier denne plotten? Hvis vi tar en log transformasjon, hvordan kan du tolke serien da?
Svar:
Det åpenbare er at husholdningers gjeld har steget kraftig.
Hvis vi tar en log transformasjon, så kan vi tolke en rett linje som konstant prosentvis endring over tid.
log_gjeld_zoo = log(gjeld_HH_zoo)
plot(log_gjeld_zoo)
Er gjeld-serien stasjonær? Hvis ikke, hvordan kan vi prøve å transformere den? (Hint: vi behøver ikke å bruke en formell test her)
Svar:
Ja, vi ser at gjeld har en positiv trend, og derfor har den ikke en konstant gjennomsnitt.
Transformere log-serien til å bli 3-måneders endring (hint: bruk funksjonen diff(serie, 3) der 4 forteller funksjonen til å ta endring over 3 måneder (kvartal)).
Hvordan kan man tolke serien?
d_gjeld_zoo = diff(log_gjeld_zoo, 3)
plot(d_gjeld_zoo)
- Er den nye transformerte serien stasjonær? Kjør en adf test med adf.test() funksjonen fra pakken tseries. Tolk resultatene.
library(tseries)## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo#adf.test krever at det ikke finnes na verdier:
d_gjeld_zoo = d_gjeld_zoo[!is.na(d_gjeld_zoo)]
adf.test(d_gjeld_zoo)##
## Augmented Dickey-Fuller Test
##
## data: d_gjeld_zoo
## Dickey-Fuller = -3.6715, Lag order = 7, p-value = 0.02637
## alternative hypothesis: stationaryEr det seriekorrelasjon i den nye serien? Bruk ACF og PACF funksjoner. Hint: Vi må først omformattere serien til vanlig ts-formatt ved å bruke funksjonen ts()
d_gjeld_ts = ts(d_gjeld_zoo)
acf(d_gjeld_ts)
pacf(d_gjeld_ts)
- Kjør en autoreggresiv modell med fire autoregressive ledd. Bruk funksjonen dynlm() fra pakken dynlm. Tolk resultatene.
Svar:
mod1 = dynlm(d_gjeld_zoo ~ L(d_gjeld_zoo, 1) + L(d_gjeld_zoo, 2) + L(d_gjeld_zoo, 3) + L(d_gjeld_zoo, 4)-1)
summary(mod1)##
## Time series regression with "zoo" data:
## Start = 1988-07-01, End = 2017-12-01
##
## Call:
## dynlm(formula = d_gjeld_zoo ~ L(d_gjeld_zoo, 1) + L(d_gjeld_zoo,
## 2) + L(d_gjeld_zoo, 3) + L(d_gjeld_zoo, 4) - 1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.063341 -0.001336 0.001245 0.003328 0.065503
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## L(d_gjeld_zoo, 1) 1.01188 0.04459 22.692 < 2e-16 ***
## L(d_gjeld_zoo, 2) -0.07197 0.06378 -1.128 0.26
## L(d_gjeld_zoo, 3) -0.54283 0.06378 -8.511 5.15e-16 ***
## L(d_gjeld_zoo, 4) 0.54801 0.04457 12.296 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.008394 on 350 degrees of freedom
## Multiple R-squared: 0.8651, Adjusted R-squared: 0.8636
## F-statistic: 561.2 on 4 and 350 DF, p-value: < 2.2e-16Vi ser en sterk autokorrelasjon på det første leddet. Vi kan tolke dette som at endring i gjeld er persistent fra måned til måned. En måned med høy endring vekst er ofte følgt med en måned med høy endring i vekst.
De andre leddene kan være vanskeligere å tolke, men siden vi tok endring fra kvartal til kvartal, er det antakeligvis noe sesong-mønster som vi ser.
Er det seriekorrelasjon igjen i residualene? Hva har det å si om resultatene?
mod1_resids = mod1$residuals
acf(ts(mod1_resids))
pacf(ts(mod1_resids))
- ACF og PACF viser forholdsvis lite seriekorrelasjon i
- Hadde det vært seriekorrelasjon, så er det mulig at standard-feilene ikke var riktig.
Oppgave II: Valutakurs og inflasjon (KPI)
I denne oppgaven skal vi laste ned data på både inflasjon (KPI) og valutakurs og sjekke forholdet mellom de.
# laster inn data på KPI
#tallene viser prosentvis endring fra samme måned året før
kpi = read.csv("http://jmaurit.github.io/norwayeconomy/data_series/kpi.csv")
head(kpi)## date KPI KPI.JAE KPIXE trimmet_snitt vektet_median
## 1 2006-01-01 1.8 0.8 1.3 1.6 1.7
## 2 2006-02-01 2.6 1.0 1.4 1.7 1.8
## 3 2006-03-01 2.3 0.8 1.4 1.5 1.6
## 4 2006-04-01 2.7 0.8 1.3 1.5 1.8
## 5 2006-05-01 2.3 0.8 1.2 1.5 1.6
## 6 2006-06-01 2.2 0.8 1.3 1.5 1.7kpi["date"] = as.Date(kpi$date, format="%Y-%m-%d")
#Vi laster inn kronekurs data
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"]
#få valutakurs-formatt
kronekurs["USD_NOK"] <- 1/kronekurs["NOK_USD"]Her setter vi sammen dataseriene med Merge
kpi_valuta = kronekurs %>% inner_join(kpi, by="date")Vi plotter KPI og KPI justert for avgiftsendring og energi priser.
ggplot(kpi_valuta) +
geom_line(aes(x=date, y=KPI, color="red")) +
geom_line(aes(x=date, y=KPI.JAE, color="blue")) +
scale_color_discrete(name="",
labels=c("KPI", "Ex. energi og avgift"))
Her ser vi på valutakurset mot dollaren.
ggplot(kpi_valuta) +
geom_line(aes(x=date, y=USD_NOK))
Endelig, kan vi formattere til Zoo
kpi_valuta_zoo = zoo(kpi_valuta[c("USD_NOK", "KPI.JAE")], order.by=kpi_valuta$date)I Lab 12 så vi at USD-NOK ikke var stasjonær. Er KPI stasjonær?
Hva kan vi transformere seriene og hvordan kan vi tolke seriene?
Hva er mest riktig: å analysere valutakurset (endogen variabel) som en funksjon av KPI (eksogen variabel). Eller å analysere KPI (endogen) som en funksjon av KPI?
Kjør en dynamisk (distribuerte lag) regresjon og tolk resultantene.
Finnes det seriekorrelasjon igjen i residualene?
library(tseries)
adf.test(kpi_valuta_zoo$KPI.JAE)##
## Augmented Dickey-Fuller Test
##
## data: kpi_valuta_zoo$KPI.JAE
## Dickey-Fuller = -2.3915, Lag order = 5, p-value = 0.4127
## alternative hypothesis: stationaryVi kan ikke forkaste null-hypotesen at serien er ikke-stasjonær.
Vi kan transformere begge seriene til log-endring
log_endring_kpi = diff(log(kpi_valuta_zoo$KPI.JAE))
log_endring_vk = diff(log(kpi_valuta_zoo$USD_NOK))plot(log_endring_vk)
adf.test(log_endring_kpi)## Warning in adf.test(log_endring_kpi): p-value smaller than printed p-value##
## Augmented Dickey-Fuller Test
##
## data: log_endring_kpi
## Dickey-Fuller = -4.7154, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationaryadf.test(log_endring_vk)## Warning in adf.test(log_endring_vk): p-value smaller than printed p-value##
## Augmented Dickey-Fuller Test
##
## data: log_endring_vk
## Dickey-Fuller = -4.5968, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationarydyn_mod1 = dynlm(log_endring_kpi ~ lag(log_endring_kpi) + log_endring_vk + lag(log_endring_vk) + lag(log_endring_vk), 2)
summary(dyn_mod1)##
## Time series regression with "numeric" data:
## Start = 1, End = 172
##
## Call:
## dynlm(formula = log_endring_kpi ~ lag(log_endring_kpi) + log_endring_vk +
## lag(log_endring_vk) + lag(log_endring_vk), data = 2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.428e-16 -3.040e-17 -1.730e-17 -5.900e-18 3.304e-15
##
## Coefficients: (1 not defined because of singularities)
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.164e-17 1.969e-17 5.910e-01 0.555
## lag(log_endring_kpi) 1.000e+00 9.367e-17 1.068e+16 <2e-16 ***
## log_endring_vk 4.538e-16 7.297e-16 6.220e-01 0.535
## lag(log_endring_vk) NA NA NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.571e-16 on 169 degrees of freedom
## Multiple R-squared: 1, Adjusted R-squared: 1
## F-statistic: 5.784e+31 on 2 and 169 DF, p-value: < 2.2e-16- En 1% endring i valutakurset fører til en -1.5 % endring i kpi.
- Valutakurset øker: billigere med import - lavere kpi.
library(lmtest)
bgtest(dyn_mod1)##
## Breusch-Godfrey test for serial correlation of order up to 1
##
## data: dyn_mod1
## LM test = 1.1386, df = 1, p-value = 0.2859Kan ikke forkaste null-hypotesen av ingen serie-korrelasjon.
Vi kan se hva som skjer hvis vi ikke har med AR(1) leddet i modellen:
dyn_mod2 = dynlm(log_endring_kpi ~ log_endring_vk + lag(log_endring_vk) + lag(log_endring_vk), 2)
summary(dyn_mod2)##
## Time series regression with "numeric" data:
## Start = 1, End = 172
##
## Call:
## dynlm(formula = log_endring_kpi ~ log_endring_vk + lag(log_endring_vk) +
## lag(log_endring_vk), data = 2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.10086 -0.09420 -0.00157 0.09914 1.05626
##
## Coefficients: (1 not defined because of singularities)
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.005396 0.016117 0.335 0.738
## log_endring_vk -0.941600 0.593100 -1.588 0.114
## lag(log_endring_vk) NA NA NA NA
##
## Residual standard error: 0.2105 on 170 degrees of freedom
## Multiple R-squared: 0.01461, Adjusted R-squared: 0.008813
## F-statistic: 2.52 on 1 and 170 DF, p-value: 0.1142bgtest(dyn_mod2)##
## Breusch-Godfrey test for serial correlation of order up to 1
##
## data: dyn_mod2
## LM test = 13.77, df = 1, p-value = 0.0002066Oppgave III: Strømpriser og vindkraft.
I denne oppgaven skal vi prøve å estimere hvordan produksjon av vindkraft påvirker strømpriser.
Vi begynner med å laste inn dataserier for strømpriser og vindkraft produksjon i norden:
#laster inn vindkraft produksjon fra danmark og sverige (i MWh)
vind_prod =read.csv("http://jmaurit.github.io/norwayeconomy/data_series/dk_se_wind_data.csv")
# "spot" prisen på strøm i ulike soner (EUR/MWh)
elspot = read.csv("http://jmaurit.github.io/norwayeconomy/data_series/elspot.csv")
#omformatterer til dato-formatt
elspot["date"] = as.Date(elspot$date, format="%Y-%m-%d")
vind_prod["date"] = as.Date(vind_prod$date, format="%Y-%m-%d")Vi begrenser vår analyse av område som heter DK1: hovedsakelig halvøyen Jylland i danmark og bruker merge til å sette sammen vår data-serier.
DK1 = merge(vind_prod[c("DK1", "date")], elspot[c("DK1", "date")], by="date")
colnames(DK1)[2:3] = c("vind_kraft", "elpris")
tail(DK1)## date vind_kraft elpris
## 67 2018-07-01 11245 5214
## 68 2018-08-01 2489 5569
## 69 2018-09-01 3017 4989
## 70 2018-10-01 32970 4748
## 71 2018-11-01 31766 5404
## 72 2018-12-01 30001 4670Kjør en vanlig linear regresjon (bruk lm()) av vind kraft på elpris, tolk resultatene.
Finnes det AR(1) seriekorrelasjon igjen i residualene? Gjør en formell test med en Breusch-Godfrey test.
lm_mod1 = lm(elpris ~ vind_kraft, data=DK1)
summary(lm_mod1)##
## Call:
## lm(formula = elpris ~ vind_kraft, data = DK1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1997.3 -487.6 -109.9 351.3 2225.5
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.429e+03 1.721e+02 19.922 <2e-16 ***
## vind_kraft -7.892e-03 5.323e-03 -1.483 0.143
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 862.1 on 70 degrees of freedom
## Multiple R-squared: 0.03044, Adjusted R-squared: 0.01659
## F-statistic: 2.198 on 1 and 70 DF, p-value: 0.1427- med en koeffisient som er -.028 sier dette at hvis vindkraft øker med 100 MWh, så blir prisen redusert med -2.8 EUR/MWh.
- Dette er signifikant.
Seriekorrelasjon? :
lm_mod1_res = lm_mod1$residuals
acf(lm_mod1_res)
* Det ser absolutt ut som det er residualer, men vi kan gjøre en formell test:
library(lmtest)
bgtest(lm_mod1, order=1)##
## Breusch-Godfrey test for serial correlation of order up to 1
##
## data: lm_mod1
## LM test = 48.69, df = 1, p-value = 2.997e-12*Med en LM-statistikk med 675 og p-verdi av nærmest 0, så kan vi forkaste vår null-hypotese (ingen seriekorrelasjon).
Nå omformattere til zoo og test om begge seriene er stasjonær
Bruk dynlm til å kjøre en distribuerte lag modell med en AR(1)-ledd og 2 lagger av vindkraft-leddet. Tolk koeffisientene.
Er det fortsatt seriekorrelasjon i residualene? Hvilket betydning har det?
Svar: Vi formaterer vår dataserier som zoo:
DK1_zoo = zoo(DK1[c("elpris", "vind_kraft")], order.by=DK1$date)plot(DK1_zoo$elpris)
plot(DK1_zoo$vind_kraft)
library(lmtest)
adf.test(DK1_zoo$elpris)##
## Augmented Dickey-Fuller Test
##
## data: DK1_zoo$elpris
## Dickey-Fuller = -0.75576, Lag order = 4, p-value = 0.9611
## alternative hypothesis: stationaryadf.test(DK1_zoo$vind_kraft)##
## Augmented Dickey-Fuller Test
##
## data: DK1_zoo$vind_kraft
## Dickey-Fuller = -4.0869, Lag order = 4, p-value = 0.01076
## alternative hypothesis: stationary- Begge ser stasjonær ut
dlm_mod1 = dynlm(elpris~vind_kraft + L(vind_kraft) + L(vind_kraft, 2) + L(elpris), data=DK1_zoo)
summary(dlm_mod1)##
## Time series regression with "zoo" data:
## Start = 2013-03-01, End = 2018-12-01
##
## Call:
## dynlm(formula = elpris ~ vind_kraft + L(vind_kraft) + L(vind_kraft,
## 2) + L(elpris), data = DK1_zoo)
##
## Residuals:
## Min 1Q Median 3Q Max
## -992.37 -260.80 -17.65 309.78 1045.60
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 723.538609 271.480544 2.665 0.0097 **
## vind_kraft -0.002429 0.002885 -0.842 0.4029
## L(vind_kraft) -0.003787 0.002934 -1.291 0.2012
## L(vind_kraft, 2) -0.001466 0.002949 -0.497 0.6207
## L(elpris) 0.839535 0.066439 12.636 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 451.6 on 65 degrees of freedom
## Multiple R-squared: 0.7469, Adjusted R-squared: 0.7313
## F-statistic: 47.96 on 4 and 65 DF, p-value: < 2.2e-16vindkraft har fortsatt en negative effect på priser samme dag.
vi estimerer en positiv koeffisient på første lagg av vind kraft. Sansynligvis er dette en effekt som kommer fra at vindkraft er autokorrelert (man burde være forsiktig med å tolke det som at vindkraft fra forrige dag påvirker dagen’s elpris.)
Vi ser også at elprisen er positivt autokorrelert: En sjokk i elprisen vil føres videre gjennom tid.
Sjekker for seriekorrelasjon i residualene:
dlm_mod1_res = dlm_mod1$residuals
acf(ts(dlm_mod1_res))
bgtest(dlm_mod1)##
## Breusch-Godfrey test for serial correlation of order up to 1
##
## data: dlm_mod1
## LM test = 6.8774, df = 1, p-value = 0.008729- Ser ut som det fortsatt er noe seriekorrelasjon (kanskje sesong-relatert.)
- Dette kan føre til at standard-feilene kan være feil.
En løsning er å kalkulere “robust” standard feil (noenganger kalt Newey West) fra pakken Sandwich.
library(sandwich)
coeftest(dlm_mod1, vcov = NeweyWest)##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 723.5386091 169.1980418 4.2763 6.359e-05 ***
## vind_kraft -0.0024291 0.0026763 -0.9076 0.3674
## L(vind_kraft) -0.0037874 0.0025904 -1.4621 0.1485
## L(vind_kraft, 2) -0.0014662 0.0020699 -0.7083 0.4813
## L(elpris) 0.8395354 0.0535572 15.6755 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- Nå er det kun vind_kraft i samme perioden som er signifikant!