Wstęp

W niniejszym projekcie wykorzystuję model Bayesian Auto-Regressive (zwany dalej BAR), aby przewidywać wartości inflacji w Polsce na najbliższe miesiące.

Korzystam ze sharmonizowanej inflacji konsumenckiej HICP dostępnej na stronie Eurostatu dla lat 1997-2023 (https://ec.europa.eu/eurostat/databrowser/view/prc_hicp_manr/default/table?lang=en). Najpierw na podstawie danych z lat 1997-2022 wyliczę prognozę dla 11 miesięcy roku 2023 i porównam ją z wartościami faktycznymi. Następnie, korzystając z tych samych paramterów wyliczę przewidywania inflacji dla kilku kolejnych miesięcy.

Kody

Kody, z których korzystam do szacowania modelu BAR w większości bazują na slajdach i kodach opracowanych przez doktora Andrzeja Kocięckiego z Uniwersytetu Warszawskiego (https://scholar.google.com/citations?user=svf8CzoAAAAJ&hl=pl&oi=ao). Kody odpowiadające za graficzną część niniejszej pracy są wyłącznie mojego autorstwa (https://www.linkedin.com/in/tymoteusz-metrak/).

Wczytywanie bibliotek

library(tidyverse) # łatwiejsza praca z danymi
library(lubridate)  # pracowanie z datami
library(ggthemes) # bardziej estetyczne wykresy
library(dygraphs) # interaktywe wykresy
library(xts) # dla klasy "xts" koniecznej do interaktywnych wykresów

Wczytywanie danych

inflation <- read.csv("prc_hicp_manr_page_linear_march_2024.csv")
Poland <- inflation[inflation$geo=="PL", 7:8]
Poland$TIME_PERIOD <- ym(Poland$TIME_PERIOD) # zmiana wyrazów na typ "data"

Pierwszy rzut oka na dane - zbuduję wykres, aby zobaczyć inflację w rzeczonym okresie. Najpierw tworzę ramkę danych z zapisanymi wartościami inflacji, datą i miejscem na oszacowania, które później obliczę.

df <- data.frame(
  c(rep(NA, length(Poland[,2]))), # wartości NA, potem włożę tu przewidywania
  c(Poland[,2]), # prawdziwe wartości
  c(Poland[,1])) # daty
names(df) <- c("Fitted", "Actual", "Date")

Następnie tworzę poglądowy wykres.

ggplot(data=df, aes(x=Date)) +
  geom_line(aes(y=Actual, colour = "Inflacja r/r"), linewidth=1) +
  theme_gdocs() + 
  theme(legend.title = element_blank())+
  ylab("Inflacja HICP (w %)") + 
  xlab("Rok") +
  labs(title="Inflacja w Polsce 1997-2023",
       subtitle = "styczeń 1997 - luty 2024")

Wyświetlę teraz wykres interkatywny.

Poland_xts <- xts(Poland[c("OBS_VALUE")], order.by = Poland$TIME_PERIOD)
dygraph(Poland_xts$OBS_VALUE,
        main = "Wykres interaktywny inflacji w Polsce 1997-2023",
        xlab="Rok",
        ylab="HICP (w %)") %>%
  dyRangeSelector(height = 60)

Model przewidujący inflację w 2023 roku

Poniższa część pracy to model, który na podstawie danych 1997-2022 przewidzi inflację HICP w 11 miesiącach 2023 roku. Następnie poównam wartości przewidywane z wartościami rzeczywistymi, aby wykazać użyteczność modelu BAR do przewidywania inflacji i ogólniej - szeregów czasowych.

Zapisywanie danych w wektorze

yy <- as.matrix(Poland[1:312,2]) # wykorzystam dane 1997-2022 i będę przewidywał na 2023 rok

T.obs <- length(yy) ## "T" -> długość powyższego wektora

Parametry

Ustalam 12 wartości opóźnionych, które odpowiadają całemu roku poprzedzającmu dany miesiąc. Horyzont prognozy ustalam na 12 miesięcy 2023 roku.

p <- 12      ## liczba opóźnień (lags)
h <- 12      ## horyzont prognozy

Technikalia I

Konieczne transformacje macierzy

y <- cbind(yy[(p+1):T.obs])

x <- matrix(NA,T.obs-p,p+1)


x[,1] <- 1 


for (i in 1:p) {
  x[,i+1] <- yy[(p+1-i):(T.obs-i)]
}

Reszty z MNK

Potrzebuję otrzymać reszty z modelu regresji liniowej szacowanej metodą najmniejszych kwadratów. Użyję potem te reszty do obliczenia prawdopodobieństwa a posteriori dla sigma^2.

beta.hat = solve(crossprod(x),crossprod(x,y))
## Ekwiwalent znanego wzoru = solve(t(x)%*%x)%*%t(x)%*%y

# Reszty
z <- y - x%*%beta.hat

Rozkład a priori

Definuję tutaj odpowiedni rozkład a priori. To w tym miejscu definiuję jaki tym modelu BAR wykorzystam. Może to być BAR(1), BAR(2), albo BAR(1) z dodatkowym opóźnieniem odpowiadającym za sezonowość (“seasonality”). To na takie ostatnie rozwiązanie się zdecyduję. Do szacowania wykorzystam wartość z ostatniego okresu i sprzed roku (gdyż inflacja podlega wahaniom sezonowym, np. z powodu wzrostu konsumpcji w okresie świątecznym czy spadku bezrobocia latem).

Parametr “shrinkage” ustalam na poziomie 0,15, chociaż może być to zmieniane w zależności od potrzeb.

beta.bar = matrix(0,p+1,1)
beta.bar[2] <- 1        ## random walk prior
beta.bar[13] <- 1        ## seasonality prior

# Parameters poniżej odpowiadające za "shrinkage"
lambda.bar <- 0.15      ## Shrinkage dla wszystkich bet poza stałą

lambda.bar.const <- 1   ## Shrinkage (a właściwie jego brak) dla stałej

Technikalia II

Macierz sigma bar

Tworzę macierz złożoną z samych zer i wartości lambda na diagonali.

Sigma.bar <- matrix(0,p+1,p+1)

# Pierwsza wartość równa lambda.bar.const (dla stałej)
Sigma.bar[1,1] <- lambda.bar.const

# Pozostałe wartości zwyczajnie na diagonali
for (i in 1:p) {
  Sigma.bar[i+1,i+1] <- lambda.bar/i
}

Sigma tilde + beta tilde

Sigma.tilde <- crossprod(x) + solve(Sigma.bar)
inv.Sigma.tilde <- solve(Sigma.tilde)
beta.tilde <- inv.Sigma.tilde%*%( crossprod(x)%*%beta.hat + solve(Sigma.bar)%*%beta.bar)

Krańcowe prawdopodobieństwo a posteriori dla sigma^2

Poniżej są parametry dla prawdopodobieństwa brzegowego sigma^2. Są to parametry z funkcji gęstości odwrotnego gamma (inverted-gamma probability density function).

alpha.sigma2 <- 0.5*(T.obs-p)
lambda.sigma2 <- 0.5*(crossprod(z) + crossprod(x%*%beta.hat) + t(beta.bar)%*%solve(Sigma.bar)%*%beta.bar - t(beta.tilde)%*%Sigma.tilde%*%beta.tilde )

Liczba iteracji

Większa liczba iteracji (powtórzeń) wydłuży proces szacowania, ale powinna zwiększyć dokładność estymacji. Ja ustalam liczbę powtórzeń algorytmu na 300 000, ale już przy 100 000 można otrzymać wiarygodne oszacowania.

Ponadto, w kroku tym tworzę miejsce na przechowywanie wartości, które zaraz oszacuję.

n.iter <- 300000

sigma2.all <- matrix(NA,n.iter,1)
sigma2.gibbs.all <- matrix(NA,n.iter,1)
beta.all <- matrix(NA,n.iter,p+1)
forecast.all <- matrix(NA,n.iter,h)

Próbkowanie - technikalia III

W próbkowaniu wykorzystam tzw. właściwość Choleskiego.

beta <- beta.bar   # inicjalizacja

for (i in 1:n.iter) {
  
  # próbkowanie z rozkładu odwrotnego gamma (rozkład brzegowy a posteriori sigma^2)
  sigma2 <- 1/rgamma(1, shape = alpha.sigma2, rate = lambda.sigma2) 
  
  # próbkowanie z wielowymiarowego rozkładu normalnego z wykorzystaniem właściwości Choleskiego 
  beta <- beta.tilde + sqrt(sigma2)*t(chol(inv.Sigma.tilde))%*%rnorm(p+1)


  # próbkowanie z rozkładu prognozy
  forecast <- matrix(NA,1,h)
  
  x.last <- x[T.obs-p,2:(p+1)] # ostatni rząd macierzy x bez 1
  
  for (j in 1:h) {
    y.mean <- cbind(1,rbind(x.last))%*%beta
    obs.for <- y.mean + sqrt(sigma2)*rnorm(1) # próbkowanie z prawdopodobieństwa dla przyszłej obserwacji
    
    forecast[j] <- obs.for
    
    x.last <- cbind(obs.for,rbind(x.last))
    x.last <- x.last[1:p]
  }
  sigma2.all[i] <- sigma2
  beta.all[i,] <- t(beta)
  forecast.all[i,] <- forecast 
}

Wspaniale! Próbkowanie trwało kilka minut i zostało zakończone. Otrzymałem już wyniki, teraz należy je odpowiednio zapisać.

W macierzy forecast.all mam wszystkie wartości próbkowania. Macierz ta ma 300 000 wierszy - po jednym dla każdego powtorzenia algorytmu próbkowania. Posiada też 12 kolumn, odpowiadających 12 okresom przewidywania inflacji. Dla każdej kolumny wyliczam medianę, dolne 5% i górne 5% (czyli 95%).

sr.quant <- matrix(NA,h,3)

for (i in 1:h) {
  sr.quant[i,] <- quantile(forecast.all[,i], probs = c(0.05,0.5,0.95)) 
}

Wartości przewidywane vs. rzeczywiste

Teraz, do stworzonej na początku tabeli z wartościami prawdziwymi dodam przewidywania (mediana) wartościami inflacji dla 12 miesięcy 2023 roku. Pozwoli mi to na wyliczenie błędu standardowego - MAPE (mean absolute percentage error) oraz przedstawienie przewidywań na wykresie.

df$Fitted <- c(rep(NA, 312), sr.quant[,2], NA, NA) # NA values and predictions

Błąd standardowy MAPE (mean absolute percentage error)

Tutaj wyliczam tzw. “Średni absolutny błąd procentowy”.

MAPE <- filter(df, year(Date)>2022 & year(Date)<2024) %>% 
  summarise(MAPE=mean(abs(Actual-Fitted)/Actual))
MAPE <- MAPE[1,1]
print(MAPE)
## [1] 0.0533613

Jak widać, zadany model produkuje niezbyt duży błąd standardowy na przestrzeni 12 miesięcy 2023 roku i można go uznać za satysfakcjonujący. Zaakceptuję zatem wcześniej zdefiniowany rozkład a priori, liczbę opóźnień (lags) oraz shrinkage. Można jednak manipulować tymi wartościami w zależności od potrzeb (i nadziei na zmniejszenie błędu prognozy).

90% przedział ufności prognozy

W tym miejscu utworzę tabelę z przedziałem ufności na poziomie 90%.

posterior.interval <- cbind.data.frame(
  sr.quant[,1],
  sr.quant[,3], 
  subset(df, year(Date)>2022 & year(Date)<2024)$Date)
names(posterior.interval) <- c("LL", "UL", "Date")

Teraz połączę przedziałów wiarygodności z medianą przewidywań.

credible_df <- left_join(df, posterior.interval, by="Date")

Wykresy

Najpierw stowrzę wykres porównujący inflację w całym badanym wykresie z prognozą na ostatnie miesiące 2023 roku.

ggplot(data=credible_df, aes(x=Date)) +
  geom_line(aes(y=Actual, colour = "Prawdziwa"), linewidth=1) +
  geom_line(aes(y=Fitted, colour = "Przewidywana"), linewidth=1, linetype=2) +
  theme_gdocs() + 
  theme(legend.title = element_blank()) + 
  ylab("") + 
  xlab("") +
  geom_vline(xintercept=as.numeric(as.Date("2022-12-01")), linetype=2) + 
  geom_ribbon(aes(ymin=LL, ymax=UL), fill="grey", alpha=0.5) +
    labs(title=paste0("Przewidywana vs. prawdziwa inflacja w Polsce"),
       subtitle="1997 - 2023",
       caption=paste0("Mean Absolute Percentage Error = ", round(100*MAPE,2), "%","\nAutor: Tymoteusz Mętrak"),
       x="Data",
       y="Punkty procentowe (%)")
## Warning: Removed 314 rows containing missing values (`geom_line()`).

Teraz stworzę wykres wyłącznie dla miesięcy, dla których mogę porównać wartości przewidywane z prawdziwymi.

credible_df_new <- credible_df[313:324,]
ggplot(data=credible_df_new, aes(x=Date)) +
  geom_line(aes(y=Actual, colour = "Prawdziwa"), linewidth=1) +
  geom_line(aes(y=Fitted, colour = "Przewidywana"), linewidth=1, linetype=2) +
  theme_gdocs() + 
  theme(legend.title = element_blank()) + 
  ylab("") + 
  xlab("") +
  geom_ribbon(aes(ymin=LL, ymax=UL), fill="grey", alpha=0.5) +
  labs(title=paste0("Przewidywana vs. prawdziwa inflacja w Polsce"),
       subtitle="Styczeń - grudzień 2023",
       caption=paste0("Mean Absolute Percentage Error = ", round(100*MAPE,2), "%","\nAutor: Tymoteusz Mętrak"),
       x="Data",
       y="Punkty procentowe (%)")

Wykres interaktywny z wartościami przewidywanymi

Zbuduję teraz wykres z “suwakiem”.

credible_df_xts <- xts(credible_df[,-3], order.by = credible_df$Date)
credible_df_xts[312,1]=credible_df_xts[312,2]
credible_df_xts[312,3]=credible_df_xts[312,2]
credible_df_xts[312,4]=credible_df_xts[312,2]
dygraph(credible_df_xts, 
        main = "Wykres interaktywny inflacji w Polsce 1997-2023",
        xlab="Rok",
        ylab="HICP (w %)") %>%
  dySeries("Actual", label = "Prawdziwa") %>%
  dySeries(c("LL", "Fitted", "UL"), label = "Przewidywana") %>%
  dyRangeSelector(height = 70)

Przewidywanie inflacji na przyszłość

W tej sekcji dokonam przewidywań inflacji na kolejny rok. Skorzystam ze sprawdzonych powyżej parametrów, które zagwarantowały niski błąd standardowy.

Zapisywanie danych w wektorze

yy2 <- as.matrix(Poland[,2]) # wykorzystam wszystkie dane 1997-luty 2024

T.obs2 <- length(yy2)

Parametry

Ustalam 12 wartości opóźnionych, które odpowiadają całemu roku poprzedzającmu dany miesiąc. Horyzont prognozy ustalam na 12 miesięcy.

p2 <- 12      # liczba opoźnień
h2 <- 12      # horyznot prognozy

Technikalia IV

Potrzebne transformacje

y2 <- cbind(yy2[(p2+1):T.obs2])

x2 <- matrix(NA,T.obs2-p2,p2+1)


x2[,1] <- 1  


for (i in 1:p2) {
  x2[,i+1] <- yy2[(p2+1-i):(T.obs2-i)]
}

Reszty z MNK

beta.hat2 = solve(crossprod(x2),crossprod(x2,y2))

# Reszty
z2 <- y2 - x2%*%beta.hat2

Rozkład a priori

Definuję tutaj identyczny rozkład a priori jak wcześniej. Parametr “shrinkage” ustalam na poziomie 0,15, chociaż jak wspomniałem, może być to zmieniane w zależności od potrzeb.

beta.bar2 = matrix(0,p2+1,1)
beta.bar2[2] <- 1        ## random walk prior
beta.bar2[13] <- 1        ## seasonality prior

# "Shrinkage"
lambda.bar2 <- 0.15 

lambda.bar.const2 <- 1

Technikalia V

Macierz sigma bar

Sigma.bar2 <- matrix(0,p2+1,p2+1)

Sigma.bar2[1,1] <- lambda.bar.const2

for (i in 1:p2) {
  Sigma.bar2[i+1,i+1] <- lambda.bar2/i
}

Sigma tilde + beta tilde

Sigma.tilde2 <- crossprod(x2) + solve(Sigma.bar2)
inv.Sigma.tilde2 <- solve(Sigma.tilde2)
beta.tilde2 <- inv.Sigma.tilde2%*%( crossprod(x2)%*%beta.hat2 + solve(Sigma.bar2)%*%beta.bar2)

Rozkład brzegowy a posteriori sigma^2

alpha.sigma22 <- 0.5*(T.obs2-p2)
lambda.sigma22 <- 0.5*(crossprod(z2) + crossprod(x2%*%beta.hat2) + t(beta.bar2)%*%solve(Sigma.bar2)%*%beta.bar2 - t(beta.tilde2)%*%Sigma.tilde2%*%beta.tilde2 )

Liczba iteracji

Ustalam liczbę powtórzeń algorytmu tak jak wcześniej na 300 000.

Ponadto, w kroku tym tworzę miejsce na przechowywanie wartości, które zaraz oszacuję.

n.iter2 <- 300000

sigma2.all2 <- matrix(NA,n.iter2,1)
beta.all2 <- matrix(NA,n.iter2,p2+1)
forecast.all2 <- matrix(NA,n.iter2,h2)

Próbkowanie

beta2 <- beta.bar2

for (i in 1:n.iter2) {
  
  sigma22 <- 1/rgamma(1, shape = alpha.sigma22, rate = lambda.sigma22)
  
  beta2 <- beta.tilde2 + sqrt(sigma22)*t(chol(inv.Sigma.tilde2))%*%rnorm(p+1)

  
  forecast2 <- matrix(NA,1,h2)
  
  x.last2 <- x2[T.obs2-p2,2:(p2+1)]
  
  for (j in 1:h2) {
    y.mean2 <- cbind(1,rbind(x.last2))%*%beta2
    obs.for2 <- y.mean2 + sqrt(sigma22)*rnorm(1)
    
    forecast2[j] <- obs.for2
    
    x.last2 <- cbind(obs.for2,rbind(x.last2))
    x.last2 <- x.last2[1:p2]
  }
  sigma2.all2[i] <- sigma22
  beta.all2[i,] <- t(beta2)
  forecast.all2[i,] <- forecast2
}

Próbkowanie zostało zakończone. Otrzymałem już wyniki, teraz należy je odpowiednio zapisać. Wyliczam medianę, dolne 5% i górne 5% (czyli 95%).

sr.quant2 <- matrix(NA,h2,3)

for (i in 1:h2) {
  sr.quant2[i,] <- quantile(forecast.all2[,i], probs = c(0.05,0.5,0.95)) 
}

Przygotuję ramkę danych z wartościami przewidywanymi

dates <- seq(as.Date("2024/03/01"), as.Date("2025/02/01"), by = "month")
df2 <- data.frame(sr.quant2[,2],dates)
names(df2) <- c("Fitted", "Date")

90% przedział ufności

W tym miejscu utworzę tabelę z przedziałem ufności na poziomie 90%.

posterior.interval2 <- cbind.data.frame(
  sr.quant2[,1],
  sr.quant2[,3], 
  df2$Date)
names(posterior.interval2) <- c("LL", "UL", "Date")
credible_df2 <- left_join(df2, posterior.interval2, by="Date")

Wykresy

Najpierw wykres z inflacją przewidywaną na przyszłość.

ggplot(data=credible_df2, aes(x=Date)) +
  geom_line(aes(y=Fitted, colour = "Przewidywana"), linewidth=1, linetype=2) +
  theme_gdocs() + 
  theme(legend.title = element_blank()) + 
  ylab("") + 
  xlab("") +
  geom_ribbon(aes(ymin=LL, ymax=UL), fill="grey", alpha=0.5) +
    labs(title="Przewidywana inflacja w Polsce",
       subtitle="marzec 2024 - luty 2025",
       caption="Autor: Tymoteusz Mętrak",
       x="Data",
       y="Punkty procentowe (%)")

Stworzę kolejne wykresy. Do tego potrzebuję zbudować odpowiednie tabele.

dates2 <- seq(as.Date("1997/01/01"), as.Date("2025/02/01"), by = "month")
df3 <- data.frame(
  c(rep(NA, 326), sr.quant2[,2]),
  c(Poland$OBS_VALUE, rep(NA, 12)),
  dates2)
names(df3) <- c("Fitted", "Actual","Date")
df3 <- left_join(df3, posterior.interval2, by="Date")

Teraz mogę stworzyć wykres z wartościami prawdziwymi oraz przewidywanymi.

ggplot(data=df3, aes(x=Date)) +
  geom_line(aes(y=Actual, colour = "Prawdziwa"), linewidth=1) +
  geom_line(aes(y=Fitted, colour = "Przewidywana"), linewidth=1, linetype=2) +
  theme_gdocs() + 
  theme(legend.title = element_blank()) + 
  ylab("") + 
  xlab("") +
  geom_vline(xintercept=as.numeric(as.Date("2024-02-01")), linetype=2) + 
  geom_ribbon(aes(ymin=LL, ymax=UL), fill="grey", alpha=0.5) +
    labs(title=paste0("Inflacja w Polsce i jej przewidywanie"),
       subtitle="styczeń 1997 - luty 2025",
       caption="Autor: Tymoteusz Mętrak",
       x="Data",
       y="Punkty procentowe (%)")
## Warning: Removed 12 rows containing missing values (`geom_line()`).
## Warning: Removed 326 rows containing missing values (`geom_line()`).

Finalnie, stworzę jeszcze ostatni wykres - wykres interaktywny.

df3_xts <- xts(df3[,-3], order.by = df3$Date)
df3_xts[326,1] <- df3_xts[326,2]
df3_xts[326,3] <- df3_xts[326,2]
df3_xts[326,4] <- df3_xts[326,2]
dygraph(df3_xts, 
        main = "Wykres interaktywny inflacji w Polsce 1997-2025",
        xlab="Rok",
        ylab="HICP (w %)") %>%
  dySeries("Actual", label = "Prawdziwa") %>%
  dySeries(c("LL", "Fitted", "UL"), label = "Przewidywana") %>%
  dyRangeSelector(height = 70)