Analizy predyktyne z paczkami bsvars i bsvarSIGNs dla R

\[ \]

cechy paczek bsvars i bsvarSIGNs

strukturalne modele VAR

modelowanie rozkładu i zmienności

analizy predyktywne

skrypty dla analiz predyktywnych

materiały

\[ \]

slajdy jako strona internetowa

repozytorium na GitHub dla reprodukcji wyników

skrypty dla reprodukcji wyników

\[ \]

bsvars.org officjalna strona

paczka bsvars na stronach CRAN

paczka bsvarSIGNs na stronach CRAN

cechy paczek bsvars i bsvarSIGNs

cechy paczek bsvars i bsvarSIGNs

bsvars.org

cechy paczek bsvars i bsvarSIGNs

bsvars.org

\[ \]

  • bayesowska estymacja modeli strukturalnych VAR
  • koherentna struktura kodu, skryptów i objektów
  • świetna szybkość obliczeniowa
  • najnowsze metody ekonometryczne i numeryczne
  • napisane w C++ dzięki paczkom Rcpp i RcppArmadillo
  • analiza danych w R

cechy paczek bsvars i bsvarSIGNs

bsvars.org

  • ładowanie paczki i danych
library(bsvars)
data(us_fiscal_lsuw)
  • łatwa inicjalizacja modelu
spec = specify_bsvar$new(us_fiscal_lsuw)
  • prosta estymacja
burn = estimate(spec, S = 1000)
post = estimate(burn, S = 10000)

  • ładowanie paczki i danych
library(bsvarSIGNs)
data(optimism)
  • łatwa inicjalizacja modelu
spec = specify_bsvarSIGN$new(optimism)
  • prosta estymacja
post = estimate(spec, S = 10000)

cechy paczek bsvars i bsvarSIGNs

bsvars.org

  • analizy strukturalne
irfs = compute_impulse_responses(post , horizon = 12)
fevd = compute_variance_decompositions(post, horizon = 12)
hds  = compute_historical_decompositions(post)
ss   = compute_structural_shocks(post)
csds = compute_conditional_sd(post)
sddr = verify_identification(post)

  • analizy strukturalne
irfs = compute_impulse_responses(post , horizon = 12)
fevd = compute_variance_decompositions(post, horizon = 12)
hds  = compute_historical_decompositions(post)
ss   = compute_structural_shocks(post)
csds = compute_conditional_sd(post)

cechy paczek bsvars i bsvarSIGNs

bsvars.org

  • analizy predyktywne
fvs  = compute_fitted_values(post)
fore = forecast(post, horizon = 12)
  • wykresy i podsumowania
plot(irfs)
summary(irfs)

  • analizy predyktywne
fvs  = compute_fitted_values(post)
fore = forecast(post, horizon = 12)
  • wykresy i podsumowania
plot(irfs)
summary(irfs)

cechy paczek bsvars i bsvarSIGNs

bsvars.org

  • skrypty z przekierowaniem
library(bsvars)
data(us_fiscal_lsuw)

us_fiscal_lsuw |> 
  specify_bsvar$new() |> 
  estimate(S = 1000) |> 
  estimate(S = 10000) -> post

post |> compute_impulse_responses(horizon = 12) |> plot()
post |> compute_variance_decompositions(horizon = 12) |> plot()
post |> compute_historical_decompositions() |> plot()
post |> compute_structural_shocks() |> plot()
post |> compute_conditional_sd() |> plot()
post |> forecast(horizon = 12) |> plot()
post |> verify_identification() |> summary()

  • skrypty z przekierowaniem
library(bsvarSIGNs)
data(optimism)

optimism |> 
  specify_bsvarSIGN$new() |> 
  estimate(S = 10000) -> post

post |> compute_impulse_responses(horizon = 12) |> plot()
post |> compute_variance_decompositions(horizon = 12) |> plot()
post |> compute_historical_decompositions() |> plot()
post |> compute_structural_shocks() |> plot()
post |> compute_conditional_sd() |> plot()
post |> forecast(horizon = 12) |> plot()

cechy paczek bsvars i bsvarSIGNs

bsvars.org

  • monitorowanie postępu

  • monitorowanie postępu

strukturalne modele VAR

strukturalne modele VAR

  • podstawowe dla modelowania efektów polityki ekonomicznej
  • analiza dynamicznych efektów przyczynowych dobrze izolowanej przyczyny
  • stosunkowo proste w pracy z danymi i dostarczają empirycznych dowodów na propagację szoków przez gospodarki i rynki
  • dostarczają empirycznych faktów do uwzględnienia w modelach teoretyczne
  • szeroko stosowane w: polityce pieniężnej i fiskalnej, rynku finansowym, …
  • rozszerzalne: wiele wariantów specyfikacji
    • nieliniowość
    • heteroskedastyczność
    • zmienne parametry w czasie
    • modelowanie hierarchiczne bayesowskie
  • zaproponowane przez Sims (1980)

bsvars.org

strukturalne modele VAR

bsvars.org

model.

\[\begin{align} \text{równanie VAR: }&& y_t &= \mathbf{A}_1 y_{t-1} + \dots + \mathbf{A}_p y_{t-p} + \mathbf{A}_d x_{t} + \epsilon_t\\[1ex] \text{równanie strukturalne: }&& \mathbf{B}\epsilon_t &= u_t\\[1ex] \text{strukturalne szoki: }&& u_t |Y_{t-1} &\sim N_N\left(\mathbf{0}_N,\text{diag}\left(\boldsymbol\sigma_t^2\right)\right) \end{align}\]

notacja.

  • \(y_t\) - wektor \(N\) zmiennych na okres \(t\)
  • \(\mathbf{A}_i\) i \(\mathbf{B}\) - \(N\times N\) macierze parametrów autoregresyjnych i strukturalnych
  • \(\epsilon_t\) i \(u_t\) - wektory \(N\) błędów statystycznych i szoków strukturalnych
  • \(\boldsymbol\sigma_t^2\) - wektor \(N\) wariancji szoków strukturalnych

SVAR: hierarchiczne rozkłady a priori

bsvars.org

  • normalny-uogólniony normalny rozkład a priori dla \(\mathbf{A}\) i \(\mathbf{B}\)
  • wielopoziomowa estymacja wariancji a priori
  • rozkład a priori z Minnesoty dla niestacjonarnych szeregów czasowych
  • bardziej precyzyjne estymacja i prognozowanie

  • rozkład a priori normalny i odwrócony Wisharta dla \(\mathbf{A}\) i \(\mathbf{\Sigma} = (\mathbf{B}'\mathbf{B})^{-1}\)
  • estymacja wariancji a priori
  • rozkład a priori z Minnesoty dla niestacjonarnych szeregów czasowych
  • bardziej precyzyjne estymacja i prognozowanie

SVAR: modelowanie zmienności

bsvars.org

  • homoskedastyczność \(\boldsymbol\sigma_{n.t}^2 = 1\)

  • zmienność stochastyczna

  • stacjonarny proces Markowa dla zmienności

  • nieparametryczny proces Markowa dla zmienności

  • rozkłady szoków

    • normalny
    • skończona mieszanka rozkładów normalnych
    • nieparametryczna mieszanka rozkładów normalnych
    • rozkład t-Studenta

  • homoskedastyczność
  • normalny rozkład szoków

SVAR: identyfikacja

bsvars.org

  • restrykcje zerowe
  • heteroskedastyczność
  • nienormalne rozkłady szoków

  • restrykcje znaków
  • restrykcje zerowe
  • restrykcje narracyjne

identyfikacja modeli strukturalnych

analizy predyktywne

cel prognozowania ekonomicznego

bsvars.org

\(\left.\right.\)

… to efektywne użycie danych do dostarczenia statystycznej charakterystyki nieznanych przyszłych wartości interesujących dla nas wielkości.

\(\left.\right.\)

Pełna statystyczna charakterystyka nieznanych przyszłych wartości zmiennych losowych jest dana przez ich gęstość predykcyjną

\(\left.\right.\)

Proste charakterystyki gęstośc predykcyjnej są zwykle używane w procesach decyzyjnych.

\(\left.\right.\)

Charakterystyki gęstości predykcyjnej są również komunikowane do szerszych odbiorców.

gęstość predykcyjna w modelach SVAR

bsvars.org

model SVAR.

\[\begin{align*} y_t &= \mathbf{A}_1 y_{t-1} + \dots + \mathbf{A}_p y_{t-p} + \mathbf{A}_d x_{t} + \epsilon_t\\[2ex] \epsilon_t|Y_{t-1} &\sim iidN_N\left(\mathbf{0}_N,\mathbf\Sigma_t\right)\\[2ex] \mathbf\Sigma_t &= \mathbf{B}^{-1} \text{diag}\left(\boldsymbol\sigma^2_t\right)\mathbf{B}^{-1\prime} \end{align*}\]

warunkowa gęstość predykcyjna na jeden okres.

\[\begin{align*} y_{t+h}|Y_{t+h-1},\mathbf{A},\mathbf\Sigma_{t+h} &\sim N_N\left(\mathbf{A}_1 y_{t+h-1} + \dots + \mathbf{A}_p y_{t+h-p} + \mathbf{A}_d x_{t+h},\mathbf\Sigma_{t+h}\right) \end{align*}\]

gęstość predykcyjna

bsvars.org

\(\left.\right.\)

Prognozowanie bayesowskie uwzględnia niepewność co do estymacji parametrów, wycałkowując ją z gęstości predykcyjnej.

\[\begin{align*} &\\ p(y_{T+1}|Y,X) &= \int p(y_{T+1}|Y_{T},\mathbf{A},\mathbf\Sigma_{T+1}) p(\mathbf{A},\mathbf\Sigma_{T+1}|Y,X) d(\mathbf{A},\mathbf\Sigma_{T+1})\\ & \end{align*}\]

  • \(p(y_{T+1}|Y,X)\) - gęstość predykcyjna
  • \(p(y_{T+1}|Y_{t},\mathbf{A},\mathbf\Sigma_{T+1})\) - warunkowa gęstość predykcyjna na jeden okres
  • \(p(\mathbf{A},\mathbf\Sigma_{T+1}|Y,X)\) - brzegowy rozkład a posteriori parameterów
  • uwzględnia predykcję zmienności \(\boldsymbol\sigma^2_{T+1}\) dla obliczenia \(\mathbf\Sigma_{T+1}\)

próbkowanie z gęstości predykcyjnej

bsvars.org

\(\left.\right.\)

krok 1: próbkuj z rozkładu a posteriori parametrów

… aby otrzymać \(S\) ciągnięć \(\left\{ \mathbf{A}^{(s)},\mathbf\Sigma_{T+1}^{(s)} \right\}_{s=1}^{S}\)

\(\left.\right.\)

krok 2: próbkuj z gęstości predykcyjnej

Aby uzyskać ciągnięcia z \(p(y_{T+1}|Y,X)\), dla każdego z \(S\) ciągnięć z \((\mathbf{A},\mathbf\Sigma_{T+1})\) próbkuj odpowiadające wartości \(y_{T+1}\):

Próbkuj \(y_{T+1}^{(s)}\) z \[ N_N\left(\mathbf{A}_1^{(s)} y_{T} + \dots + \mathbf{A}_p^{(s)} y_{T-p+1} + \mathbf{A}_d^{(s)} x_{T+1},\mathbf\Sigma_{T+1}^{(s)}\right) \] aby otrzymać \(\left\{y_{T+1}^{(s)}\right\}_{s=1}^{S}\)

impulse response functions

bsvars.org

\(\left.\right.\)

Impulse response functions dla szoków strukturalnych obliczone dla właściwie sformułowanego empirycznie modelu SVAR są rozważane jako dynamiczne efekty przyczynowe szoków \(u_t\) na zmienne ekonomiczne \(y_{t+i}\) na \(i\) okresów po wystąpieniu szoku.

\[\begin{align*} \frac{\partial y_{n.t+i}}{\partial u_{j.t}}&=\theta_{nj.i} \end{align*}\]

  • \(\theta_{nj.i}\) - impulse response \(n\)-tej zmiennej na \(j\) szok \(i\) okresów po wystąpieniu szoku dla \(i=0,1,\dots,h\) i \(n,j=1,\dots,N\)

impulse response functions

bsvars.org

\(\left.\right.\)

…jako różnica w prognozach warunkowych.

\[\begin{align} \mathbf\Theta_{h[\cdot i]} &= \mathbb{E}\left[\mathbf{y}_{T+h}| u_{i.T} = 1,\mathbf{Y},\mathbf{X}\right] - \mathbb{E}\left[\mathbf{y}_{T+h}| \mathbf{Y},\mathbf{X}\right] \end{align}\]

…jako mechanizm propagacji szoku.

\[\begin{align} y_t &= \mathbf{B}^{-1}\mathbf{u}_t + \mathbf{A}_1 y_{t-1}\\ &= \mathbf{B}^{-1}\mathbf{u}_t + \mathbf{A}_1 \mathbf{B}^{-1}\mathbf{u}_{t-1} + \mathbf{A}_1^2 y_{t-2}\\ &\dots\\ &= \mathbf{\Theta}_0\mathbf{u}_t + \mathbf{\Theta}_1\mathbf{u}_{t-1} + \mathbf{\Theta}_2\mathbf{u}_{t-2} + \dots \end{align}\]

forecast error variance decompositions

bsvars.org

\[ \]

…są sposobem na szacunek względnej ważności szoków w wyjaśnianiu wariancji predyktywnej zmiennych w systemie.

\[ \]

macierz kowariancji predyktywnej.

\[\begin{align} \mathbf\Sigma_{T+h} &= \mathbf{\Theta}_h \text{diag}\left(\boldsymbol\sigma^2_{T+h}\right)\mathbf{\Theta}_h^{\prime} \end{align}\]

wariancja predyktywna.

Dla pierwszej zmiennej w okresie \(h\) naprzód, dla \(N=2\):

\[\begin{align} \sum_{i=0}^{h-1}\mathbf\Sigma_{T+i[11]} &= \sum_{i=0}^{h-1} \mathbf{\Theta}_{i[11]}^2 \sigma^2_{1.T+i} + \mathbf{\Theta}_{i[12]}^2 \sigma^2_{2.T+i} \end{align}\]

forecast error variance decompositions

bsvars.org

\[ \]

…są sposobem na szacunek względnej ważności szoków w wyjaśnianiu wariancji predyktywnej zmiennych w systemie.

\[ \]

forecast error variance decomposition

..dla wpływu pierwszego szoku na wariancję predyktywną pierwszej zmiennej:

\[\begin{align} FEVD_{T+h[11]} &= \sum_{i=0}^{h-1}\frac{\mathbf{\Theta}_{i[11]}^2 \sigma^2_{1.T+i}}{ \mathbf{\Theta}_{i[11]}^2 \sigma^2_{1.T+i} + \mathbf{\Theta}_{i[12]}^2 \sigma^2_{2.T+i}} \end{align}\]

modelowanie rozkładu i zmienności

niescentrowana zmienność stochastyczna

bsvars.org

\[\begin{align} &\\ \text{wariancja warunkowa:}&&\sigma_{n.t}^2 &= \exp\left\{\omega_n h_{n.t}\right\}\\ \text{w skali log:}&&h_{n.t} &= \rho_n h_{n.t-1} + v_{n.t}\\ \text{innowacje zmienności:}&&v_{n.t}&\sim N\left(0,1\right)\\ \end{align}\]

  • świetna zdolność do prognozowania
  • normalizajca \(\sigma_{n.t}^2 = 1\)
  • verify_identification() przez ocene restrykcji \(H_0:\omega_n = 0\)

scentrowana zmienność stochastyczna

bsvars.org

\[\begin{align} &\\ \text{wariancja warunkowa:}&&\sigma_{n.t}^2 &= \exp\left\{ \tilde{h}_{n.t}\right\}\\ \text{w skali log:}&&\tilde{h}_{n.t} &= \rho_n \tilde{h}_{n.t-1} + \tilde{v}_{n.t}\\ \text{innowacje zmienności:}&&\tilde{v}_{n.t}&\sim N\left(0,\omega_n^2\right)\\ \end{align}\]

  • świetna zdolność do prognozowania

zmienność stochastyczna

rozkłady brzegowe a priori.

proces Markowa dla zmienności.

bsvars.org

\[\begin{align} &\\ \text{szoki strukturalne:}&&\mathbf{u}_t\mid s_t \sim N\left( \mathbf{0}_N, \text{diag}\left(\boldsymbol{\sigma}_{s_t}^2\right) \right)\\ \text{a priori:}&& M^{-1}\left(\boldsymbol{\sigma}_{1}^2, \dots, \boldsymbol{\sigma}_{M}^2\right) \sim Dirichlet(\underline{a}\boldsymbol\imath')\\ \text{proces Markowa:}&& s_t\sim \text{Markov}(\mathbf{P},\boldsymbol\pi_0) \end{align}\]

  • modelowanie proces Markowa dla zmienności
  • zapewnia identyfikację
  • poprawa zdolności do prognozowania
  • verify_identification() przez ocenę restrykcji \(H_0:\boldsymbol{\sigma}_{1}^2, \dots, \boldsymbol{\sigma}_{M}^2 = 1\)

mieszanka rozkładów normalnych.

bsvars.org

\[\begin{align} &\\ \text{szoki strukturalne:}&&\mathbf{u}_t\mid s_t \sim N\left( \mathbf{0}_N, \text{diag}\left(\boldsymbol{\sigma}_{s_t}^2\right) \right)\\ \text{a priori:}&& M^{-1}\left(\boldsymbol{\sigma}_{1}^2, \dots, \boldsymbol{\sigma}_{M}^2\right) \sim Dirichlet(\underline{a}\boldsymbol\imath')\\ \text{proces Markowa:}&& Pr[s_t]=\boldsymbol\pi_0 \end{align}\]

  • modelowanie mieszanki rozkładów normalnych
  • zapewnia identyfikację
  • verify_identification() przez ocenę restrykcji \(H_0:\boldsymbol{\sigma}_{1}^2, \dots, \boldsymbol{\sigma}_{M}^2 = 1\)

rozkład t-studenta.

bsvars.org

\[\begin{align} &&&\\ \text{szoki strukturalne:}&&\mathbf{u}_t\mid\mathbf{x}_t &\sim t\left( \mathbf{0}_N, \mathbf{I}_N, \nu \right) \end{align}\]

  • \(\nu\) - szacowane z danych stopnie swobody
  • grube ogony zapewniają identyfikację
  • poprawa zdolności do prognozowania
  • verify_identification() przez ocenę restrykcji \(H_0:\nu \rightarrow\infty\)

skrypty dla analiz predykcyjnych

mały system polityki pieniężnej

bsvars.org

\[ \]

Dla australijskich danych kwartalnych jak w Turnip (2017)

\[\begin{align} y_t = \begin{bmatrix} \Delta rgdp_t & \pi_t & cr_t & \Delta rtwi_t \end{bmatrix}' \end{align}\]

identyfikacja przez trójkątną macierz strukturalną.

\[\begin{align} \begin{bmatrix} B_{11}&0&0&0\\ B_{21}&B_{22}&0&0\\ B_{31}&B_{32}&B_{33}&0\\ B_{41}&B_{42}&B_{43}&B_{44} \end{bmatrix} \begin{bmatrix} \Delta rgdp_t \\ \pi_t \\ cr_t \\ \Delta rtwi_t \end{bmatrix} \end{align}\]

mały system polityki pieniężnej

# Gross domestic product (GDP); Chain volume
rgdp_dwnld      = readrba::read_rba(series_id = "GGDPCVGDP")
rgdp_tmp        = xts::xts(rgdp_dwnld$value, rgdp_dwnld$date, tclass = 'yearqtr')
drgdp           = na.omit(400 * diff(log(rgdp_tmp)))
drgdp           = xts::to.quarterly(drgdp, OHLC = FALSE)

# Consumer price index; All groups; Quarterly change (in per cent)
picpi_dwnld     = readrba::read_rba(series_id = "GCPIAGSAQP")
pi              = 4 * xts::xts(picpi_dwnld$value, picpi_dwnld$date, tclass = 'yearqtr')
pi              = xts::to.quarterly(pi, OHLC = FALSE)

# Interbank Overnight Cash Rate
cr_dwnld        = readrba::read_rba(series_id = "FIRMMCRID")   # Cash Rate Target
cr_tmp          = xts::xts(cr_dwnld$value, cr_dwnld$date)
cr              = xts::to.quarterly(cr_tmp, OHLC = FALSE)

# Real Trade-Weighted Index
rtwi_dwnld      = readrba::read_rba(series_id = "FRERTWI")
rtwi_tmp        = xts::xts(rtwi_dwnld$value, rtwi_dwnld$date, tclass = 'yearqtr')
rtwi            = 100 * na.omit(diff(log(rtwi_tmp)))
drtwi            = xts::to.quarterly(rtwi, OHLC = FALSE)

y               = na.omit(merge(drgdp, pi, cr, drtwi))
plot(y, main = "mały system australijskiej polityki pieniężnej",
     legend.loc = "bottomleft", col = c("#FF00FF","#990099","#0c606d","#330033"))

bsvars.org

mały system polityki pieniężnej

inicjalizacja i estymacja modelu

bsvars.org

# estimation - lower-triangular model
############################################################
library(bsvars)
set.seed(123)

spec = specify_bsvar$new(
  as.matrix(y), 
  p = 4, 
  stationary = rep(TRUE, 4)
)
spec |>
  estimate(S = 1000) |>
  estimate(S = 5000) -> post
**************************************************|
bsvars: Bayesian Structural Vector Autoregressions|
**************************************************|
 Gibbs sampler for the SVAR model                 |
**************************************************|
 Progress of the MCMC simulation for 1000 draws
    Every draw is saved via MCMC thinning
 Press Esc to interrupt the computations
**************************************************|
**************************************************|
bsvars: Bayesian Structural Vector Autoregressions|
**************************************************|
 Gibbs sampler for the SVAR model                 |
**************************************************|
 Progress of the MCMC simulation for 5000 draws
    Every draw is saved via MCMC thinning
 Press Esc to interrupt the computations
**************************************************|

oblicz impulse responses

bsvars.org

post |> 
  compute_impulse_responses(horizon = 20) |> 
  plot()

oblicz forecast error variance decompositions

bsvars.org

post |> 
  compute_variance_decompositions(horizon = 20) |> 
  plot()

oblicz szoki strukturanle

bsvars.org

post |> 
  compute_structural_shocks() |> 
  plot()

oblicz wartości dopasowane

bsvars.org

post |> 
  compute_fitted_values() |> 
  plot()

oblicz prognozy

bsvars.org

post |> 
  forecast(horizon = 8) |> 
  plot(data_in_plot = 0.3)

oblicz warunkowane prognozy

bsvars.org

cond_fore = matrix(NA, 8,4)
cond_fore[,3] = c(4.35,4.10,3.35,2.85,2.60,2.35,2.10,2.10)

post |> 
  forecast(horizon = 8, conditional_forecast = cond_fore) |> 
  plot(data_in_plot = 0.3)