Модель 1 (bsts)
library(bsts)
library(ggplot2)
library(readxl)
ds <-data.frame(read_excel("C:/Users/user/Downloads/mn.xlsx"))
Date <-as.Date(ds$Date, format="%Y-%m")
Y <- ds$GRP
mod <-list()
ggplot(ds, aes(x=Date, y=GRP)) +
geom_line() + theme_minimal()
mod <- AddLocalLinearTrend(mod, Y)
mod <- AddAutoAr(mod, Y, lag = 3)
model <- bsts(Y, mod, timestamps = Date,
niter = 10000, ping = 0, seed = 42)
summary(model)
## $residual.sd
## [1] 1.54965
##
## $prediction.sd
## [1] 2.324903
##
## $rsquare
## [1] 0.4611872
##
## $relative.gof
## [1] 0.4377083
AcfDist(residuals(model))
plot(model)
plot(model, y = "components")
forecast <-predict(model, horizon = 13, burn = 3000)
plot(forecast)
a <-plot(forecast, plot.original = 30)
mod2 <-list()
mod2 <- AddLocalLinearTrend(mod2, Y)
mod2 <- AddAutoAr(mod2, Y, lag = 2)
model2 <- bsts(Y, mod2, timestamps = Date,
niter = 10000, ping = 0, seed = 42)
summary(model2)
## $residual.sd
## [1] 1.657531
##
## $prediction.sd
## [1] 2.353949
##
## $rsquare
## [1] 0.3835552
##
## $relative.gof
## [1] 0.4240649
AcfDist(residuals(model2))
plot(model2)
plot(model2, y = "components")
forecast2 <-predict(model2, horizon = 13, burn = 3000)
plot(forecast2)
a2 <-plot(forecast2, plot.original = 30)
mod3 <-list()
mod3 <- AddLocalLinearTrend(mod3, Y)
mod3 <- AddAutoAr(mod3, Y, lag = 4)
model3 <- bsts(Y, mod3, timestamps = Date,
niter = 10000, ping = 0, seed = 42)
summary(model3)
## $residual.sd
## [1] 1.576829
##
## $prediction.sd
## [1] 2.324212
##
## $rsquare
## [1] 0.4421214
##
## $relative.gof
## [1] 0.4381162
AcfDist(residuals(model3))
plot(model3)
plot(model3, y = "components")
forecast3 <-predict(model3, horizon = 13, burn = 3000)
plot(forecast3)
a3 <-plot(forecast3, plot.original = 30)
mod4 <-list()
mod4 <- AddLocalLinearTrend(mod4, Y)
mod4 <- AddAutoAr(mod4, Y, lag = 1)
model4 <- bsts(Y, mod4, timestamps = Date,
niter = 10000, ping = 0, seed = 42)
summary(model4)
## $residual.sd
## [1] 1.748422
##
## $prediction.sd
## [1] 2.405402
##
## $rsquare
## [1] 0.3140961
##
## $relative.gof
## [1] 0.3990478
AcfDist(residuals(model4))
plot(model4)
plot(model4, y = "components")
forecast4 <-predict(model4, horizon = 13, burn = 3000)
plot(forecast4)
a4 <-plot(forecast4, plot.original = 30)
better_model <- list("model" = model, "model2" = model2,
"model3" = model3, "model4" = model4)
CompareBstsModels(better_model)
Модель 2 (условное прогнозирование)
set.seed(42)
library(BVAR)
library(readxl)
ds <-read_excel("C:/Users/user/Downloads/vars2.xlsx")
mn <- bv_minnesota(lambda = bv_lambda(mode = 0.2, sd = 0.4, min = 0.0001, max = 5),
alpha = bv_alpha(mode = 2), var = 1e07)
soc <- bv_soc(mode = 1, sd = 1, min = 1e-04, max = 50)
sur <- bv_sur(mode = 1, sd = 1, min = 1e-04, max = 50)
priors <- bv_priors(hyper = "auto", mn = mn, soc = soc, sur = sur)
mh <- bv_metropolis(scale_hess = c(0.05, 0.0001, 0.0001),
adjust_acc = TRUE, acc_lower = 0.25, acc_upper = 0.45)
run <- bvar(ds, lags = 4, n_draw = 50000, n_burn = 25000, n_thin = 1,
priors = priors, mh = mh, verbose = TRUE)
print(run)
## Bayesian VAR consisting of 118 observations, 5 variables and 4 lags.
## Time spent calculating: 23.24 secs
## Hyperparameters: lambda, soc, sur
## Hyperparameter values after optimisation: 0.32563, 0.48213, 0.31463
## Iterations (burnt / thinning): 50000 (25000 / 1)
## Accepted draws (rate): 8721 (0.349)
plot(run)
fitted(run, type = "mean")
## Numeric array (dimensions 118, 5) with fitted values from a BVAR.
## Median values:
## GRP: 1.68, 0.24, -0.28, [...], 1.34, -0.15, -1.78
## CPI: 0.68, 1.16, 1.04, [...], 0.65, 0.51, 0.46
## K.Rate: 5.78, 7.17, 7.28, [...], 13.3, 13.53, 14.54
## Urals: 105.04, 102.47, 106.47, [...], 84.17, 79.55, 71.04
## E.Rate: 35.95, 37.02, 36.56, [...], 93.05, 94.8, 87.72
plot(residuals(run, type = "mean"), vars = c("GRP", "CPI", "K.Rate", "Urals", "E.Rate"))
#bazovii prognoz po kursu
mean_value <- 90.1
num_samples <- 12
random_numbers <- rnorm(num_samples, mean = mean_value)
print(random_numbers,round(4))
## [1] 90.25 89.80 90.03 90.88 91.85 90.74 89.71 90.26 89.94 92.26 92.24 88.73
path <-c(91.47, 89.54, 90.46, 90.73, 90.50, 89.99, 91.61, 90.01, 92.12, 90.04, 91.40, 92.39)
predict(run) <- predict(run, horizon = 12, cond_path = path, cond_var = "E.Rate")
plot(predict(run), t_back = 12, "GRP")
#konservativnii prognoz po kursu
mean_value <- 95.7
num_samples <- 12
random_numbers <- rnorm(num_samples, mean = mean_value)
print(random_numbers,round(4))
## [1] 94.22 95.68 95.47 96.23 94.32 96.03 95.75 96.69 97.12 95.95 95.52 97.50
path2 <-c(96.56, 94.81, 96.49, 95.75, 94.52, 96.95, 96.95, 95.41, 95.16, 96.67, 94.92, 96.49)
predict(run) <- predict(run, horizon = 12, cond_path = path2, cond_var = "E.Rate")
plot(predict(run), t_back = 12,"GRP")
#RMSE
print(residuals(run))
## Numeric array (dimensions 118, 5) with residual values from a BVAR.
## Median values:
## GRP: -0.71, 0.11, -1.2, [...], 0.08, 0.82, 0.01
## CPI: 0.39, -0.28, -0.03, [...], 0.04, 0.18, 0.15
## K.Rate: 1.1, -0.14, 0.2, [...], -0.17, 1.42, 0.89
## Urals: 0.31, 5.42, 0.54, [...], -2.23, -6.93, -7.24
## E.Rate: 0.29, -1.32, -1.57, [...], 4.09, -4.38, 3.1
a <-residuals(run)^2
b <-colSums(a)
c <-b/122
d <-sqrt(c)
print(d)
## GRP CPI K.Rate Urals E.Rate
## 1.6168235 0.7158089 1.1031056 6.1922390 3.8557974