Проект по АВР-2

library(readxl)
library(dplyr)
library(lubridate)
library(zoo)
library(forecast)
library(vars)
library(mFilter)
library(BVAR)

Загружаю все необходимые данные У нас четыре ряда: 1. ВРП в постоянных ценах 2021 года 2. Инфляция QoQ 3. Доходы и расходы консолидированного бюджета РФ 4. Ключевая ставка

Поработаем с каждым рядом в отдельности ВРП

gdp <- read_excel("Data.xlsx", sheet = "gdp")

gdp_ts <- ts(
  gdp[[3]],
  start     = c(gdp[1,1],gdp[1,2]),
  frequency = 4
)

stl_fit <- stl(gdp_ts, s.window = "periodic")
plot(stl_fit)

gdp_des <- seasadj(stl_fit)

print(gdp_des)

             Qtr1         Qtr2         Qtr3         Qtr4
2011 2.937479e+13 2.969972e+13 2.940598e+13 2.910340e+13
2012 3.088412e+13 3.108897e+13 3.045937e+13 2.988311e+13
2013 3.120058e+13 3.160572e+13 3.093087e+13 3.072556e+13
2014 3.118719e+13 3.177528e+13 3.137630e+13 3.104034e+13
2015 3.075555e+13 3.082699e+13 3.082985e+13 3.049334e+13
2016 3.069810e+13 3.092254e+13 3.084840e+13 3.067472e+13
2017 3.106218e+13 3.162138e+13 3.168248e+13 3.102719e+13
2018 3.178392e+13 3.245500e+13 3.253722e+13 3.213606e+13
2019 3.218489e+13 3.288095e+13 3.345900e+13 3.322094e+13
2020 3.262880e+13 3.055422e+13 3.231977e+13 3.274692e+13
2021 3.281218e+13 3.367148e+13 3.370271e+13 3.454109e+13
2022 3.382183e+13 3.246435e+13 3.268951e+13 3.381724e+13
2023 3.354734e+13 3.413706e+13 3.475400e+13 3.577584e+13
2024 3.518834e+13 3.555851e+13 3.592957e+13 3.754236e+13

plot(gdp_des, main = "GDP(ts)", ylab = "GDP", xlab = "Time")

Теперь займемся рядом по бюджетным показателям

budget <- read_excel("Data.xlsx",sheet = "budget")

budget_q <- budget %>%
  arrange(year, quarter) %>%        
  group_by(year) %>%                 
  mutate(
    across(
      fed_revenue:cons_expend,      
      ~ . - lag(., default = 0),      
      .names = "{.col}_q"       
    )
  ) %>%
  ungroup()



fed_revenue_q = ts(budget_q$fed_revenue_q, start = c(2011,1), frequency = 4)
fed_expend_q  = ts(budget_q$fed_expend_q, start = c(2011,1), frequency = 4)
cons_revenue_q = ts(budget_q$cons_revenue_q, start = c(2011,1), frequency = 4)
cons_expend_q  = ts(budget_q$cons_expend_q, start = c(2011,1), frequency = 4)
  

stl_fit <- stl(cons_revenue_q, s.window = "periodic")
plot(stl_fit)

cons_revenue_des <- seasadj(stl_fit)

stl_fit <- stl(cons_expend_q, s.window = "periodic")
plot(stl_fit)

cons_expend_des <- seasadj(stl_fit)

Ключевая ставка

key_rate <- read_excel("Data.xlsx", sheet = "key_rate")

key_rate_q <- key_rate %>%
  mutate(
    year  = year(date),
    quarter = quarter(date),
    yearqtr = paste0(year, " Q", quarter)
  ) %>%
  group_by(yearqtr) %>%
  summarise(
    key_rate_q = mean(key_rate, na.rm = TRUE),
    .groups = "drop"
  )

# Создаём zoo-ряд и конвертируем в ts
key_rate_ts <- ts(key_rate_q[[2]],
                  frequency = 4,
                  start = c(2013,3) 
                  )

plot(key_rate_ts)

Инфляция QoQ

cpi <- read_excel("Data.xlsx", sheet = "cpi")

cpi_ts <- ts(
  cpi[[3]],
  start     = c(cpi[1,1],cpi[1,2]),
  frequency = 4
)
plot(cpi_ts, main = "CPI(ts)", ylab = "CPI", xlab = "Time")

stl_fit <- stl(cpi_ts, s.window = "periodic")
plot(stl_fit)

cpi_des <- seasadj(stl_fit)

plot(cpi_des,
     main = "Квартальный CPI без сезонности (STL)",
     ylab = "cpi", xlab = "Время")

# Стандартизируем масштабы
gdp     <- window(gdp_des,  start = c(2013, 3), end = c(2023, 4))/10^12
infl    <- window(cpi_des,  start = c(2013, 3), end = c(2023, 4))-100
key_rate<- window(key_rate_ts,  start = c(2013, 3), end = c(2023, 4)) 
cons_re <- window(cons_revenue_des,  start = c(2013, 3), end = c(2023, 4))/10^12
cons_ex <- window(cons_expend_des,  start = c(2013, 3), end = c(2023, 4))/10^12
balance_budget<- (-(cons_ex - cons_re)/gdp)*100

##Проведем тест на едичные корни для GDP

Единичные корни -> поиск коинтеграции + модель коррекции ошибок (ECM)

Стационарность -> векторная регрессия

ADF

gdp_adf = ur.df(gdp, type = 'trend', selectlags = 'AIC')
summary(gdp_adf)

############################################### 
# Augmented Dickey-Fuller Test Unit Root Test # 
############################################### 

Test regression trend 


Call:
lm(formula = z.diff ~ z.lag.1 + 1 + tt + z.diff.lag)

Residuals:
    Min      1Q  Median      3Q     Max 
-2.0727 -0.2883 -0.0109  0.3773  1.1263 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept) 17.24341    5.24455   3.288  0.00226 **
z.lag.1     -0.57133    0.17323  -3.298  0.00220 **
tt           0.05733    0.01724   3.324  0.00205 **
z.diff.lag   0.25562    0.17308   1.477  0.14840   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.6209 on 36 degrees of freedom
Multiple R-squared:  0.2463,    Adjusted R-squared:  0.1835 
F-statistic: 3.922 on 3 and 36 DF,  p-value: 0.01605


Value of test-statistic is: -3.2981 4.4737 5.8687 

Critical values for test statistics: 
      1pct  5pct 10pct
tau3 -4.15 -3.50 -3.18
phi2  7.02  5.13  4.31
phi3  9.31  6.73  5.61

Гипотеза об единичном корне на уровне 10% отвергается (ADF), но не 5%. Необходимо првоести теcn KPSS Гипотеза о стационарности ряда отвергается (KPSS). Ряд нестационарен

KPSS:

gdp_kpss <- ur.kpss(gdp, type="mu", lags="short")
summary(gdp_kpss)

####################### 
# KPSS Unit Root Test # 
####################### 

Test is of type: mu with 3 lags. 

Value of test-statistic is: 1.0309 

Critical value for a significance level of: 
                10pct  5pct 2.5pct  1pct
critical values 0.347 0.463  0.574 0.739

##Проведем тест на едичные корни для Inflation

infl_adf = ur.df(infl, type = 'drift', selectlags = 'AIC')
summary(infl_adf)

############################################### 
# Augmented Dickey-Fuller Test Unit Root Test # 
############################################### 

Test regression drift 


Call:
lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.7767 -0.6838 -0.2490  0.4300  4.4484 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept) 59.63801   18.67002   3.194  0.00286 **
z.lag.1     -0.58601    0.18354  -3.193  0.00287 **
z.diff.lag  -0.03715    0.16510  -0.225  0.82319   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.448 on 37 degrees of freedom
Multiple R-squared:  0.3016,    Adjusted R-squared:  0.2639 
F-statistic:  7.99 on 2 and 37 DF,  p-value: 0.001305


Value of test-statistic is: -3.1929 5.1062 

Critical values for test statistics: 
      1pct  5pct 10pct
tau2 -3.58 -2.93 -2.60
phi1  7.06  4.86  3.94

infl_kpss <- ur.kpss(infl, type="mu", lags="short")
summary(infl_kpss)

####################### 
# KPSS Unit Root Test # 
####################### 

Test is of type: mu with 3 lags. 

Value of test-statistic is: 0.1461 

Critical value for a significance level of: 
                10pct  5pct 2.5pct  1pct
critical values 0.347 0.463  0.574 0.739

Гипотеза об единичном корне на уровне 10% отвергается (ADF), но не 5%. Необходимо првоести теcn KPSS Гипотеза о стационарности ряда отвергается (KPSS). Ряд нестационарен

##Проведем тест на едичные корни для key_rate

key_rate_adf = ur.df(key_rate, type = 'none', selectlags = 'AIC')
summary(key_rate_adf)

############################################### 
# Augmented Dickey-Fuller Test Unit Root Test # 
############################################### 

Test regression none 


Call:
lm(formula = z.diff ~ z.lag.1 - 1 + z.diff.lag)

Residuals:
    Min      1Q  Median      3Q     Max 
-5.1750 -0.4280  0.0566  0.5397  8.1310 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
z.lag.1    -0.004416   0.038938  -0.113    0.910
z.diff.lag  0.189568   0.170515   1.112    0.273

Residual standard error: 2.163 on 38 degrees of freedom
Multiple R-squared:  0.03169,   Adjusted R-squared:  -0.01928 
F-statistic: 0.6218 on 2 and 38 DF,  p-value: 0.5424


Value of test-statistic is: -0.1134 

Critical values for test statistics: 
      1pct  5pct 10pct
tau1 -2.62 -1.95 -1.61

key_rate_kpss <- ur.kpss(key_rate, type="mu", lags="short")
summary(key_rate_kpss)

####################### 
# KPSS Unit Root Test # 
####################### 

Test is of type: mu with 3 lags. 

Value of test-statistic is: 0.1447 

Critical value for a significance level of: 
                10pct  5pct 2.5pct  1pct
critical values 0.347 0.463  0.574 0.739

Ряд точно нестационарный

##Проведем тест на едичные корни для Balance_Budget

autoplot(balance_budget)

budget_adf = ur.df(balance_budget, type = 'none', selectlags = 'AIC')
summary(budget_adf)

############################################### 
# Augmented Dickey-Fuller Test Unit Root Test # 
############################################### 

Test regression none 


Call:
lm(formula = z.diff ~ z.lag.1 - 1 + z.diff.lag)

Residuals:
    Min      1Q  Median      3Q     Max 
-6.2143 -0.8763 -0.3831  0.7775  5.1640 

Coefficients:
           Estimate Std. Error t value Pr(>|t|)   
z.lag.1     -0.4086     0.1371  -2.980    0.005 **
z.diff.lag   0.1567     0.1642   0.954    0.346   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.344 on 38 degrees of freedom
Multiple R-squared:  0.1911,    Adjusted R-squared:  0.1486 
F-statistic:  4.49 on 2 and 38 DF,  p-value: 0.01776


Value of test-statistic is: -2.9802 

Critical values for test statistics: 
      1pct  5pct 10pct
tau1 -2.62 -1.95 -1.61

budget_kpss <- ur.kpss(balance_budget, type="mu", lags="short")
summary(budget_kpss)

####################### 
# KPSS Unit Root Test # 
####################### 

Test is of type: mu with 3 lags. 

Value of test-statistic is: 0.0941 

Critical value for a significance level of: 
                10pct  5pct 2.5pct  1pct
critical values 0.347 0.463  0.574 0.739

По ADF ряд стационарен, по KPSS - нет.

# Разбиваем область на 2×2
par(mfrow = c(2, 2),  
    mar   = c(4, 4, 2, 1))  
# 1-й график
plot(gdp,
     main = "ВВП SA",
     xlab = "", ylab = "трлн. руб.")

# 2-й
plot(infl,
     main = "ИПЦ SA",
     xlab = "", ylab = "%")

# 3-й
plot(balance_budget,
     main = "Сальдо бюджета к ВВП",
     xlab = "Время", ylab = "%")

# 4-й
plot(key_rate,
     main = "Ключевая ставка ЦБ РФ",
     xlab = "Время", ylab = "%")

# Вернуть настройки по умолчанию (опционально)
par(mfrow = c(1,1))

###Проведем ряды ВВП, Инфляция и ключевая ставка к первым разностям

gdp_fd <- diff(log(gdp))
autoplot(gdp_fd)

key_rate_fd <-diff(key_rate)
autoplot(key_rate_fd)

infl_fd <-diff(infl)
autoplot(cons_ex_fd)

balance_budget_fd <-diff(balance_budget)
autoplot(balance_budget_fd)

#cons_re_fd <-diff(cons_re)
#autoplot(cons_re_fd)

#cons_ex_fd <-diff(cons_ex)
#autoplot(cons_ex_fd)

#data <- cbind(gdp_fd , key_rate_fd, cons_ex_fd)
data <-cbind(gdp_fd , key_rate_fd, infl_fd, balance_budget)
colnames(data) <- c("gdp_fd", "key_rate_fd", "infl_fd", "balance_budget" )

data <- window(data, start = c (2014, 1))

info_p <- VARselect(data, lag.max = 3, type = "const")
info_p$selection

AIC(n)  HQ(n)  SC(n) FPE(n) 
     1      1      1      1 

model.est <- VAR(data, p = 1, type = "const", season = NULL, exog = NULL)
summary(model.est)

VAR Estimation Results:
========================= 
Endogenous variables: gdp_fd, key_rate_fd, infl_fd, balance_budget 
Deterministic variables: const 
Sample size: 39 
Log Likelihood: -120.493 
Roots of the characteristic polynomial:
0.5944 0.5766 0.3515 0.3515
Call:
VAR(y = data, p = 1, type = "const", exogen = NULL)


Estimation results for equation gdp_fd: 
======================================= 
gdp_fd = gdp_fd.l1 + key_rate_fd.l1 + infl_fd.l1 + balance_budget.l1 + const 

                    Estimate Std. Error t value Pr(>|t|)
gdp_fd.l1         -0.0444662  0.1638993  -0.271    0.788
key_rate_fd.l1    -0.0016682  0.0020483  -0.814    0.421
infl_fd.l1        -0.0005996  0.0023463  -0.256    0.800
balance_budget.l1 -0.0017825  0.0013115  -1.359    0.183
const              0.0022592  0.0036742   0.615    0.543


Residual standard error: 0.02109 on 34 degrees of freedom
Multiple R-Squared: 0.1354, Adjusted R-squared: 0.0337 
F-statistic: 1.331 on 4 and 34 DF,  p-value: 0.2783 


Estimation results for equation key_rate_fd: 
============================================ 
key_rate_fd = gdp_fd.l1 + key_rate_fd.l1 + infl_fd.l1 + balance_budget.l1 + const 

                  Estimate Std. Error t value Pr(>|t|)  
gdp_fd.l1         36.57037   16.41898   2.227   0.0327 *
key_rate_fd.l1     0.22596    0.20519   1.101   0.2785  
infl_fd.l1        -0.06522    0.23505  -0.277   0.7831  
balance_budget.l1 -0.02511    0.13138  -0.191   0.8495  
const              0.06879    0.36807   0.187   0.8528  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


Residual standard error: 2.113 on 34 degrees of freedom
Multiple R-Squared: 0.1622, Adjusted R-squared: 0.06365 
F-statistic: 1.646 on 4 and 34 DF,  p-value: 0.1854 


Estimation results for equation infl_fd: 
======================================== 
infl_fd = gdp_fd.l1 + key_rate_fd.l1 + infl_fd.l1 + balance_budget.l1 + const 

                  Estimate Std. Error t value Pr(>|t|)    
gdp_fd.l1         16.75408   10.06806   1.664 0.105287    
key_rate_fd.l1     0.53509    0.12582   4.253 0.000156 ***
infl_fd.l1        -0.62926    0.14413  -4.366 0.000112 ***
balance_budget.l1 -0.01867    0.08056  -0.232 0.818131    
const             -0.08138    0.22570  -0.361 0.720668    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


Residual standard error: 1.295 on 34 degrees of freedom
Multiple R-Squared: 0.4819, Adjusted R-squared: 0.4209 
F-statistic: 7.905 on 4 and 34 DF,  p-value: 0.0001285 


Estimation results for equation balance_budget: 
=============================================== 
balance_budget = gdp_fd.l1 + key_rate_fd.l1 + infl_fd.l1 + balance_budget.l1 + const 

                  Estimate Std. Error t value Pr(>|t|)    
gdp_fd.l1         48.57968   17.09025   2.843  0.00752 ** 
key_rate_fd.l1    -0.11853    0.21358  -0.555  0.58256    
infl_fd.l1        -0.01273    0.24466  -0.052  0.95881    
balance_budget.l1  0.64323    0.13675   4.704 4.15e-05 ***
const             -0.53003    0.38312  -1.383  0.17554    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


Residual standard error: 2.199 on 34 degrees of freedom
Multiple R-Squared: 0.491,  Adjusted R-squared: 0.4312 
F-statistic: 8.201 on 4 and 34 DF,  p-value: 9.646e-05 



Covariance matrix of residuals:
                   gdp_fd key_rate_fd   infl_fd balance_budget
gdp_fd          0.0004447    0.005082 -0.002001        0.01152
key_rate_fd     0.0050816    4.462783  1.298407        1.73015
infl_fd        -0.0020013    1.298407  1.678050        0.58050
balance_budget  0.0115217    1.730155  0.580497        4.83515

Correlation matrix of residuals:
                 gdp_fd key_rate_fd  infl_fd balance_budget
gdp_fd          1.00000      0.1141 -0.07326         0.2485
key_rate_fd     0.11407      1.0000  0.47447         0.3725
infl_fd        -0.07326      0.4745  1.00000         0.2038
balance_budget  0.24847      0.3725  0.20379         1.0000

dum <- rep(0, dim(data)[1])
dum[25:40]<-1
plot(dum)

model.est <- VAR(data, p = 1, type = "const", season = NULL, exog = dum)
summary(model.est)

###Сделать тесты на то, что оценки OLS адекватные

model.serial <- serial.test(model.est, lags.pt = 3, type = "PT.asymptotic")
model.serial

model.arch <- arch.test(model.est, lags.multi = 3, multivariate.only = TRUE)
model.arch

    ARCH (multivariate)

data:  Residuals of VAR object model.est
Chi-squared = 326.44, df = 300, p-value = 0.1409

Гетероскедастичности также не наблюдается

model.norm <- normality.test(model.est, multivariate.only = TRUE)
model.norm

$JB

    JB-Test (multivariate)

data:  Residuals of VAR object model.est
Chi-squared = 51.056, df = 8, p-value = 2.559e-08


$Skewness

    Skewness only (multivariate)

data:  Residuals of VAR object model.est
Chi-squared = 17.551, df = 4, p-value = 0.00151


$Kurtosis

    Kurtosis only (multivariate)

data:  Residuals of VAR object model.est
Chi-squared = 33.506, df = 4, p-value = 9.409e-07

А вот нормальность остатков модели не подтвердилась