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

library(readxl)

Warning: пакет ‘readxl’ был собран под R версии 4.1.3

library(dplyr)

Warning: пакет ‘dplyr’ был собран под R версии 4.1.3

Присоединяю пакет: ‘dplyr’

Следующие объекты скрыты от ‘package:stats’:

    filter, lag

Следующие объекты скрыты от ‘package:base’:

    intersect, setdiff, setequal, union

library(dplyr)
library(lubridate)

Warning: пакет ‘lubridate’ был собран под R версии 4.1.3

Присоединяю пакет: ‘lubridate’

Следующие объекты скрыты от ‘package:base’:

    date, intersect, setdiff, union

library(zoo)

Присоединяю пакет: ‘zoo’

Следующие объекты скрыты от ‘package:base’:

    as.Date, as.Date.numeric

library(forecast)

Warning: пакет ‘forecast’ был собран под R версии 4.1.3
Registered S3 method overwritten by 'quantmod':
  method            from
  as.zoo.data.frame zoo 

library(vars)

Warning: пакет ‘vars’ был собран под R версии 4.1.3
Загрузка требуемого пакета: MASS
Warning: пакет ‘MASS’ был собран под R версии 4.1.3

Присоединяю пакет: ‘MASS’

Следующий объект скрыт от ‘package:dplyr’:

    select

Загрузка требуемого пакета: strucchange
Warning: пакет ‘strucchange’ был собран под R версии 4.1.3
Загрузка требуемого пакета: sandwich
Warning: пакет ‘sandwich’ был собран под R версии 4.1.2
Загрузка требуемого пакета: urca
Загрузка требуемого пакета: lmtest

library(mFilter)

Warning: пакет ‘mFilter’ был собран под R версии 4.1.3

    ‘mFilter’ version: 0.1-5

    ‘mFilter’ is a package for time series filtering

    See ‘library(help="mFilter")’ for details

    Author: Mehmet Balcilar, mbalcilar@yahoo.com

library(BVAR)

Присоединяю пакет: ‘BVAR’

Следующие объекты скрыты от ‘package:vars’:

    fevd, irf

Загружаю все необходимые данные У нас четыре ряда: 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 29374790 29699724 29405983 29103405
2012 30884122 31088974 30459373 29883111
2013 31200582 31605718 30930875 30725560
2014 31187188 31775282 31376305 31040339
2015 30755552 30826995 30829847 30493341
2016 30698104 30922544 30848395 30674723
2017 31062182 31621384 31682481 31027194
2018 31783917 32455003 32537220 32136064
2019 32184890 32880952 33459004 33220945
2020 32628803 30554215 32319770 32746924
2021 32812182 33671479 33702711 34541093
2022 33821834 32464353 32689511 33817238
2023 33547342 34137063 34753998 35775839
2024 35188344 35558511 35929569 37542363

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


idx <- with(budget_q, as.yearqtr(paste(year, quarter), "%Y %q"))

budget_ts <- zoo(
  cbind(
    fed_revenue_q = budget_q$fed_revenue_q,
    fed_expend_q  = budget_q$fed_expend_q,
    cons_revenue_q = budget_q$cons_revenue_q,
    cons_expend_q  = budget_q$cons_expend_q
  ),
  order.by = idx
)

plot(budget_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))
infl    <- window(cpi_des,  start = c(2013, 3), end = c(2023, 4))
cons_re <- window(as.ts(budget_ts[,3]),  start = c(2013, 3), end = c(2023, 4)) 
cons_ex <- window(as.ts(budget_ts[,4]),  start = c(2013, 3), end = c(2023, 4)) 
key_rate<- window(key_rate_ts,  start = c(2013, 3), end = c(2023, 4)) 
balance_budget <- (cons_re - cons_ex)/gdp

##Проведем тест на едичные корни для 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 
-2072694  -288320   -10902   377264  1126294 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)   
(Intercept)  1.724e+07  5.245e+06   3.288  0.00226 **
z.lag.1     -5.713e-01  1.732e-01  -3.298  0.00220 **
tt           5.733e+04  1.724e+04   3.324  0.00205 **
z.diff.lag   2.556e-01  1.731e-01   1.477  0.14840   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 620900 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)

Гипотеза об единичном корне на уровне 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 
-1.037e-04 -4.415e-05  3.766e-06  2.520e-05  7.344e-05 

Coefficients:
           Estimate Std. Error t value Pr(>|t|)    
z.lag.1     -0.8226     0.2228  -3.693 0.000694 ***
z.diff.lag  -0.1226     0.1661  -0.738 0.464963    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 4.408e-05 on 38 degrees of freedom
Multiple R-squared:  0.4638,    Adjusted R-squared:  0.4356 
F-statistic: 16.44 on 2 and 38 DF,  p-value: 7.187e-06


Value of test-statistic is: -3.6928 

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.1123 

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

Ряд стационарен по результатам обоих тестов.

# Разбиваем область на 2×2
par(mfrow = c(2, 2),    # 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))