Lab 8: Prognose med regresjon

I denne labben skal vi bruke regresjon til å skape en prognose. Denne labben er delvis basert på k.5 i Forecasting: Principles and Practice og det kan være nyttig å lese kapittelet i forkant.

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(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

I fpp2 pakken kommer det en datasett som heter uschange, som er kvartalvis prosent endringer av makroøkonomisk data fra USA. Datasettet er allerede formatert som en tidsrekke.

str(uschange)

##  Time-Series [1:187, 1:5] from 1970 to 2016: 0.616 0.46 0.877 -0.274 1.897 ...
##  - attr(*, "dimnames")=List of 2
##   ..$ : NULL
##   ..$ : chr [1:5] "Consumption" "Income" "Production" "Savings" ...

fpp2 inkluderer også en funksjon som er bygget på ggplot2, og gjør det lettere å plotte tidsrekke data. Funkjsonen heter autoplot.

autoplot(uschange[,c("Consumption","Income")]) +
  ylab("% change") + xlab("Year")

Vi kan begynne å spøre hvordan konsum er relatert til inntekt. Vi kan kjøre en enkel regresjon. fpp2 inkluderer en funksjon tslm som er litt lettere å bruke med tidsrekke data, men her fungerer det akkurat det samme som lm.

tslm(Consumption ~ Income, data=uschange)

## 
## Call:
## tslm(formula = Consumption ~ Income, data = uschange)
## 
## Coefficients:
## (Intercept)       Income  
##      0.5451       0.2806

lm(Consumption ~ Income, data=uschange)

## 
## Call:
## lm(formula = Consumption ~ Income, data = uschange)
## 
## Coefficients:
## (Intercept)       Income  
##      0.5451       0.2806

Vi kan også se på dette visuelt

uschange %>%
 as.tibble() %>%
  ggplot(aes(x=Income, y=Consumption)) +
    ylab("Consumption (quarterly % change)") +
    xlab("Income (quarterly % change)") +
    geom_point() +
    geom_smooth(method="lm")

## Warning: `as.tibble()` is deprecated as of tibble 2.0.0.
## Please use `as_tibble()` instead.
## The signature and semantics have changed, see `?as_tibble`.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.

## `geom_smooth()` using formula 'y ~ x'

#as.tibble() kommandoen gjør bare uschange til data.frame formattet "tibble" som er standard med tidyverse. 

Vi kan skape en enkel prognosemodell basert på denne regresjonen

Nå skal vi bruke pakken forecast, du må kanskje bruke install.packages(“forecast”)

library(forecast)

Vi kaller regresjonsmodellen for fit.cons

fit.cons = tslm(Consumption ~ Income, data = uschange)

Vi bruker funksjonen forecast() fra forecast pakken - du kan lese mer om den her

Før vi kan lage vår prognose, så må vi definere noen parametere:

• Vi må definere hvor lang tid i forkant vi ønsker å skape en prognose (h)

h=4

• Vi må gi info om hva vi forventer vår forklarende variabler kommer til å være (inntekt) i framtiden.

newdata = data.frame(Income = rep(mean(uschange[,"Income"]), h))
newdata

##      Income
## 1 0.7176268
## 2 0.7176268
## 3 0.7176268
## 4 0.7176268

• Og da er vi klar til å kjøre vår forecast-modell

fcast.avg = forecast(fit.cons,
  newdata, h)
fcast.avg

##         Point Forecast       Lo 80    Hi 80    Lo 95    Hi 95
## 2016 Q4      0.7464708 -0.03063903 1.523581 -0.44557 1.938512
## 2017 Q1      0.7464708 -0.03063903 1.523581 -0.44557 1.938512
## 2017 Q2      0.7464708 -0.03063903 1.523581 -0.44557 1.938512
## 2017 Q3      0.7464708 -0.03063903 1.523581 -0.44557 1.938512

Her ser vi prognose-verdiene for de neste 4 kvartalene. Vi har vår gjennomsnitts-prognose, pluss usikkerhets-intervaller. Du kan lese mer om hvordan intervallene er kalkulert her og her.

Vi kan plotte vår forecast ved å bruke autoplot, plus autolayer, der vi spesifiserer hvilken forecast modell vi har brukt.

autoplot(uschange[, "Consumption"]) +
  ylab("% endring in US Konsum") +
  autolayer(fcast.avg, PI = TRUE)

La oss i steden si at vi mener inntekt kommer til å vokse fortere de neste 10 årene.

h=10
avg_inc = mean(uschange[,"Income"])

increase = rep(1.05, h)
cum_increase = cumprod(increase)
Income = avg_inc*cum_increase

newdata = data.frame(Income = Income)
newdata

##       Income
## 1  0.7535082
## 2  0.7911836
## 3  0.8307427
## 4  0.8722799
## 5  0.9158939
## 6  0.9616886
## 7  1.0097730
## 8  1.0602617
## 9  1.1132747
## 10 1.1689385

• Hvor mye øker vekstraten i dette eksempelet?
• Forklar hva kommandoene incrase = rep(1.03, 10) og cumprod(increase) gjør?

#vi lager en ny data.frame som representerer våre forklarende verdier i framtiden. Her setter vi % i inntekt til sine gjennomsnittsverdier

fcast.up <- forecast(fit.cons,
  newdata = newdata)
autoplot(uschange[, "Consumption"]) +
  ylab("% change in US consumption") +
  autolayer(fcast.up, series = "Vekst Scenario",
    PI = TRUE) +
  guides(colour = guide_legend(title = "Scenario"))

For en prognose modell, så ønsker vi kanskje flere forklarende variabler. Vi har også tilgang til industriell produksjon, sparing og arbeidsledighet.

Vi kan få en visualisering om hvordan alle disse seriene henger sammen ved å bruke pakken GGally (du må kjøre install.packages(“GGally”), hvis du ikke har brukt pakken før) som også bygger på ggplot. Vi bruker funksjonen ggpairs, som gir korrelasjoner mellom de ulike variablene.

library(GGally)

## Registered S3 method overwritten by 'GGally':
##   method from   
##   +.gg   ggplot2

## 
## Attaching package: 'GGally'

## The following object is masked from 'package:fma':
## 
##     pigs

uschange %>%
  as.tibble() %>%
  ggpairs()

mean(uschange[,"Production"])

## [1] 0.50806

Vi kan kjøre en regresjonsmodell med tslm

fit.consBest <- tslm(
  Consumption ~ Income + Production + Unemployment + Savings,
  data=uschange)
summary(fit.consBest)

## 
## Call:
## tslm(formula = Consumption ~ Income + Production + Unemployment + 
##     Savings, data = uschange)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.88296 -0.17638 -0.03679  0.15251  1.20553 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   0.26729    0.03721   7.184 1.68e-11 ***
## Income        0.71449    0.04219  16.934  < 2e-16 ***
## Production    0.04589    0.02588   1.773   0.0778 .  
## Unemployment -0.20477    0.10550  -1.941   0.0538 .  
## Savings      -0.04527    0.00278 -16.287  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3286 on 182 degrees of freedom
## Multiple R-squared:  0.754,  Adjusted R-squared:  0.7486 
## F-statistic: 139.5 on 4 and 182 DF,  p-value: < 2.2e-16

Nå kan vi skape en prognose for de neste fire årene.

Vi må først fastsette verdier for de tre eksogene variablene. Vi skal generere to senarioer: en “up” og en “down” scenario:

h = 4

newdata = data.frame(
    Income = c(1, 1, 1, 1),
    Production = c(1,1,1,1),
    Savings = c(0.5, 0.5, 0.5, 0.5),
    Unemployment = c(0, 0, 0, 0))

fcast.up = forecast(fit.consBest, newdata = newdata)

newdata = data.frame(
    Income = rep(-1, h),
    Production = rep(0,h),
    Savings = rep(-0.5, h),
    Unemployment = rep(0, h))

fcast.down = forecast(fit.consBest, newdata = newdata)

autoplot(uschange[, 1]) +
  ylab("% change in US consumption") +
  autolayer(fcast.up, PI = TRUE, series = "increase") +
  autolayer(fcast.down, PI = TRUE, series = "decrease") +
  guides(colour = guide_legend(title = "Scenario"))

Oppgaver

• Vanlige regresjonsmodeller kan være nyttige for å skape senario-prognoser. Men hvis målet er å skape det Forecasting boken kaller “ex-ante” modeller, så er slike modeller begrenset. Forklar hvorfor.

Under kan du laste inn data om norsk BNP og andre tall fra det norsk nasjonalregnskapet.

reelt_BNP = read_csv2("http://jmaurit.github.io/makro/labs/data/reelt_norsk_BNP_NOK2015.csv")

## Using ',' as decimal and '.' as grouping mark. Use read_delim() for more control.

## Parsed with column specification:
## cols(
##   Year = col_double(),
##   konsum = col_double(),
##   off_konsum = col_double(),
##   investering = col_double(),
##   eksport = col_double(),
##   import = col_double(),
##   BNP = col_double(),
##   BNP_per_pers = col_double(),
##   folkemengde = col_double()
## )

reelt_BNP_ts = ts(reelt_BNP, freq=1, start=1830)

• Undersøk forholdet mellom investering, konsum, bnp og import over tid

• Kjør en regresjon der import er avhengig-variablen og bnp, konsum og investering er eksogene variabler. Tolk og forklar resultatene

• Skap en prognose model for import 10 år framover basert på regresjonsmodellen.

• Kan du se noen problemer med regresjonsmodellen?

reelt_BNP["import_per_BNP"] = reelt_BNP["import"]/reelt_BNP["BNP"]

reelt_BNP_l = reelt_BNP %>% select(Year, konsum, investering, import, BNP) %>% gather(konsum, investering, import, BNP, key="variabel", value="value")

ggplot(reelt_BNP_l, aes(x=Year, y=log(value))) + 
  geom_line() +
  facet_wrap(~variabel, ncol=2, scale="free_y")

ggplot(reelt_BNP, aes(x=Year, y=import_per_BNP))+
  geom_line()

ggplot(reelt_BNP, aes(x=log(investering), y=log(import), color=Year)) +
  geom_point() +
  geom_abline(slope=1, intercept=0) +
  scale_colour_gradientn(colours=terrain.colors(4))

ggplot(reelt_BNP, aes(x=log(konsum), y=log(import), color=Year)) +
  geom_point() +
  geom_abline(slope=1, intercept=0) +
  scale_colour_gradientn(colours=terrain.colors(4))

fit.cons = tslm(
  import ~ BNP + konsum + investering,
  data=reelt_BNP_ts)
summary(fit.cons)

## 
## Call:
## tslm(formula = import ~ BNP + konsum + investering, data = reelt_BNP_ts)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -53298 -17941   8141  15711  48711 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -2.630e+04  2.257e+03 -11.655  < 2e-16 ***
## BNP         -7.904e-02  2.109e-02  -3.748 0.000239 ***
## konsum       6.958e-01  6.311e-02  11.025  < 2e-16 ***
## investering  2.959e-01  8.196e-02   3.610 0.000395 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 21720 on 184 degrees of freedom
## Multiple R-squared:  0.9879, Adjusted R-squared:  0.9877 
## F-statistic:  5000 on 3 and 184 DF,  p-value: < 2.2e-16

#Beregn gjennomsnitts vekst i konsum og investering

vekst_konsum = mean(diff(reelt_BNP_ts[,"konsum"])/reelt_BNP_ts[c(1: length(reelt_BNP$konsum)-1), "konsum"])
vekst_investering = mean(diff(reelt_BNP_ts[,"investering"])/reelt_BNP_ts[c(1: length(reelt_BNP$konsum)-1), "investering"])
vekst_BNP = mean(diff(reelt_BNP_ts[,"BNP"])/reelt_BNP_ts[c(1: length(reelt_BNP$konsum)-1), "BNP"])

h = 10

#siste datapunkt
siste_konsum = reelt_BNP$konsum[length(reelt_BNP$konsum)]
siste_investering = reelt_BNP$investering[length(reelt_BNP$investering)]
siste_BNP = reelt_BNP$BNP[length(reelt_BNP$BNP)]


vekst_konsum_s = rep(vekst_konsum, h)
cum_konsum = cumprod(vekst_konsum_s)
konsum = siste_konsum*cum_increase

vekst_investering_s = rep(vekst_investering, h)
cum_investering = cumprod(vekst_investering_s)
investering = siste_investering*cum_increase

vekst_BNP_s = rep(vekst_BNP, h)
cum_BNP = cumprod(vekst_BNP_s)
BNP = as.double(siste_BNP)*cum_increase


#gjennomsnittsvekstrater

newdata = data.frame(
    BNP = BNP,
    konsum = konsum,
    investering = investering)

fcast = forecast(fit.cons, newdata = newdata)

autoplot(reelt_BNP_ts[, "import"]) +
  ylab("Import") +
  autolayer(fcast, series = "BAU Scenario",
    PI = TRUE) +
  scale_y_log10()