VECM_modeling_and_IRF.Rmd

---
title: "TS"
author: "Krasnukhina"
date: "2022-12-09"
output: html_document
---

```{r setup, include=FALSE}
library(zoo)
library(xts)
library(knitr)
library(PerformanceAnalytics)
library(ggplot2)
library(lmtest)
library(forecast)
library(dplyr)
library(tseries)
library(urca)
library("vars")
library(AER)
```

```{r}
data(Canada)
C = as.xts(Canada)
M = C[, c(1,3,4)]
```

#1 Pretest your variables and their first differences for stationarity.

Check your first variable = e for stationarity. Explain (write variable name, test name, p-value, conclusion)

```{r}
autoplot(M[,1])
autoplot(M[,2])
autoplot(M[,3])
#возможна десятилетняя сезонность в u
```

```{r}
# Apply Phillips-Perron test
pp.test(M[,1])
#e
#p-value = 0.7735, следовательно нулевая гипотеза о нестационарности ряда не отклоняется на любом адекватном уровне значимости. Ряд нестационарен. Нужно проверить первые разности.
```
```{r}
# Apply KPSS test. Interpret results.
kpss.test(M[,1])
```
Результаты теста KPSS: уровень значимости 1%, p-value меньше 0.01 (p-value smaller than printed p-value) <0.01 значит отклоняем Н0 о стационарности => ряд не стационарен (e)

```{r}
# Apply Phillips-Perron test
pp.test(M[,2])
#rw
#p-value = 0.8468 = большой, следовательно нулевая гипотеза о нестационарности ряда не отклоняется на любом уровне значимости. Ряд нестационарен. Нужно проверить первые разности.
```
```{r}
# Apply KPSS test. Interpret results.
kpss.test(M[,2])
```
Результаты теста KPSS: уровень значимости 1%, p-value меньше 0.01 (p-value smaller than printed p-value) <0.01 значит отклоняем Н0 => ряд не стационарен (rw)

```{r}
# Apply Phillips-Perron test
pp.test(M[,3])
#u
#p-value = 0.6957 = большой, следовательно нулевая гипотеза о нестационарности ряда не отклоняется на любом уровне значимости. Ряд нестационарен. Нужно проверить первые разности.
```
```{r}
# Apply KPSS test. Interpret results.
kpss.test(M[,3])
```
Результаты теста KPSS: уровень значимости 1%, p-value больше 0.1 значит не отклоняется Н0 => ряд стационарен (u)


Что-то интересное:)
Кажется, что есть и сезонность десятилетняя в u, и тренд, и дисперсия не константа
Реализуем необходимые преобразования. Сначала пробуем просто первую разность
(спойлер, пробовала и логдифф тоже, но в целом и первой разности хватает, тесты дают стационарность )

```{r}
# First difference
diff_M = diff.xts(M, 1, 1)
```

```{r}
#remove first NA value
diff_M = diff_M[-1, ]
```

```{r}
#remove outliers
t1 = tsoutliers(as.ts(diff_M[,1]),iterate=5, lambda = NULL)
diff_M[t1$index] = t1$replacements

t2 = tsoutliers(as.ts(diff_M[,2]),iterate=5, lambda = NULL)
diff_M[t2$index] = t2$replacements

t3 = tsoutliers(as.ts(diff_M[,3]),iterate=5, lambda = NULL)
diff_M[t3$index] = t3$replacements
```

```{r}
autoplot(diff_M)
autoplot(diff_M[,1])
autoplot(diff_M[,2])
autoplot(diff_M[,3])
```


```{r}
# Apply Phillips-Perron test
pp.test(diff_M[,1])
#e
```
Результаты теста PP: уровень значимости 5%, p-value = 0.01233 <0.05 значит отклоняем Н0 => ряд стационарен 
```{r}
# Apply KPSS test. Interpret results.
kpss.test(diff_M[,1])
```
Результаты теста KPSS: уровень значимости 1%, p-value=0.1>0.01 значит не можем отклонить Н0 => ряд стационарен (u)

Общий вывод о стационарности переменной diff e:
Исходя из 2 проведенных тестов, переменная является стационарной.

```{r}
# Apply Phillips-Perron test
pp.test(diff_M[,2])
#rw
```
Результаты теста PP: уровень значимости 1%, p-value меньше 0.01 (p-value smaller than printed p-value) <0.01 значит отклоняем Н0 => ряд стационарен 

```{r}
# Apply KPSS test. Interpret results.
kpss.test(diff_M[,2])
```
Результаты теста KPSS: уровень значимости 1%, p-value<0.01 значит отклоняетс Н0 => ряд не стационарен (rw). Результаты двух тестов расходятся, проведем дополнительно третий.

```{r}
# Apply simple augmented Dickey-Fuller. Interpret results. 
x<-ur.df(diff_M[,2], type = "none", selectlags = "AIC")
summary(x)
```
Результаты теста ADF: тип теста - "none", уровень значимости 1%, tau1=-2,6, статистика -2,9136 => статистика<квантиля, Н0 отклоняем, ряд стационарен.

Общий вывод о стационарности переменной diff rw:
Исходя из 2 проведенных тестов, переменная является стационарной.

```{r}
# Apply Phillips-Perron test
pp.test(diff_M[,3])
#u
```

Переменная - diff U. Результаты теста PP: уровень значимости 1%, p-value меньше 0.01 (p-value smaller than printed p-value) <0.01 значит отклоняем Н0 о нестационарности => ряд стационарен 

```{r}
# Apply KPSS test. Interpret results.
kpss.test(diff_M[,3])
```

Результаты теста KPSS: уровень значимости 1%, p-value=0.1>0.01 значит не можем отклонить Н0 о стационарности => ряд стационарен (u)

Общий вывод о стационарности переменной diff u:
Исходя из 2 проведенных тестов, переменная является стационарной.






#2 Choose optimal number of lags for VECM model
Строим вар ин левелс, определяем p и для векма возьмем на единицу меньше
```{r}
VARselect(M, lag.max = 20, type = "none", season = NULL, exogen = NULL)
```
20 лагов я в гробу видала, надо сто про брать 2 или 3, по правилу чем меньше, тем лучше для этих моделей, возьмем 2, тогда for VECM p=1

```{r}
m1 = VAR(M, p = 2, type = "const")
summary(m1)
```
```{r}
Acf(residuals(m1))
```
```{r}
plot(irf(m1,impulse = "e",ci=0.95, boot =TRUE, ortho=FALSE))
plot(irf(m1,impulse = "rw",ci=0.95, boot =TRUE, ortho=FALSE))
plot(irf(m1,impulse = "U",ci=0.95, boot =TRUE, ortho=FALSE))
```

```{r}
#cointegration test
library(urca)
jotest=ca.jo(M, type="eigen", K=2,
ecdet="none", spec="longrun")
summary(jotest)
```
K	- The lag order of the series (levels) in the VAR.
Интерпретация: 
Идем снизу с нулевого. Это не пивэлью, посмотри, что ужесточение границ идет с увеличением, если значение больше критических = > нулевая отклоняется

Переменная - e. Результаты теста Phillips-Perron, p-value = 0.7735, следовательно нулевая гипотеза о нестационарности ряда не отклоняется на любом адекватном уровне значимости. Ряд нестационарен. Результаты теста KPSS: уровень значимости 1%, p-value меньше 0.01 (p-value smaller than printed p-value) <0.01 значит отклоняем Н0 о стационарности => ряд нестационарен 

Переменная - diff U. Результаты теста PP: уровень значимости 1%, p-value меньше 0.01 (p-value smaller than printed p-value) <0.01 значит отклоняем Н0 о нестационарности => ряд стационарен. Результаты теста KPSS: уровень значимости 1%, p-value=0.1>0.01 значит не можем отклонить Н0 о стационарности => ряд стационарен. 

В рамках анализа с помощью PP, ADF, KPSS тестов было выявлено, что переменные до трансформации нестационарны, а после расчета первых разностей - стационарны. Кроме того, тесты Йохансена на коинтеграцию показали наличие коинтеграции. Все вышеупомянутое позволяет сделать вывод о том, что модель VECM может и должна быть построена. 

Для определения ранга коинтеграции были проведены trace и eigen тесты Йохансена. Результаты eigen теста: Для r=0 тестовая статистика = 43.73, а критическое значение на 1% уровне значимости = 25.75, следовательно нулевая гипотеза о том, что r = 0, отклоняется на 1% уровне значимости = > r = 0+1.  Для r = 1 нулевая гипотеза о том, что r = 1, не отклоняется на 1 процентном уровне, так как тестовая статистика = 16.46, а критическое значение = 19.19 следовательно на этом уровне значимости r = 1.

```{r}
#cointegration test
library(urca)
jotest=ca.jo(M, type="trace", K=2,
ecdet="none", spec="longrun")
summary(jotest)
```

здесь нулевая о том что r = то значение, которое ты проверяешь, альтернативная, что строго больше

интерпретация аналогичная: на 5 процентной уроыне значимости нулевая для r=1 не отклоняется = > у нас r=1

Результаты trace теста: Для r=0 тестовая статистика = 60.35, а критическое значение на 1% уровне значимости = 37.22, следовательно нулевая гипотеза о том, что r = 0, отклоняется на 1% уровне значимости = > r строго больше 0.  Для r = 1 нулевая гипотеза о том, что r = 1, не отклоняется на 5% уровне, так как тестовая статистика = 16.61, а критическое значение = 17.95, следовательно на 5% уровне значимости r = 1.

Can VECM model be applied?
Коинтеграция есть по двум тестам. Переменные не стационарны и стационарны с первой разностью. Поэтому векм релевантен.

Choose number of cointegration relation, explain (write fully test name, which H0 is rejected?) тут как раз интерпретация тестов и выбор r. 

Для выбора количества лагов было использовано правило: для построения модели VECM необходимо использовать на один лаг меньше(из-за первых разностей), чем оптимальное количество лагов для модели VAR на исходных переменных, которое определяется с помощью информационных критериев и функции VARselect. Для модели VAR на исходных переменных оптимальным значением количества лагов было выбрано - 2. По информационным критериям лучшими были варианты спецификаций с 20, 2 и 3 лагами, однако также к сведению было принято количество наблюдений в данных и потенциальное количество оцениваемых коэффициентов при построении модели, что определило невозможность оценивания модели с более чем 2 лагами на исследуемых данных, поэтому мы остановились на минимальном возможном значении = 2 лага для VAR и 1 - для VECM соответственно.

```{r}
#VECM
library("tsDyn")
VECmodel_const = VECM(M, lag = 1, r = 1, include = "const",
estim = "ML")
summary(VECmodel_const)
summary(VECmodel_const)$coefMat
```
из конитеграционного вектора получается долгосрочное  уравнение, можно менять столбцы местами, чтобы каждую переменную по очереди получать с единичкой и уравнение записывать. Дальше идут коэффициенты, где есть ECT это альфы, на листочке прописана интерпретация. Они интерпретируются так: Если меньше нуля и значимы => после краткосрочного шока возвращаемся к равновесию в долгосроке. Если больше нуля и значимы => не возвращаемся к равновесию после шоков = плохо:(

Долгосросчное равновесие может быть представлено в виде уравнения: e = -0.16 * rw +  3.94 * U. Параметр альфа для уравнения rw  = -0.1309*** (значим на 1% уровне) < 0, следовательно ряд сходящийся = в долгосрочном периоде после появления шока происходит возврат к равновесию. Параметр альфа для уравнения U  = -0.0312*** (значим на 1% уровне) < 0, следовательно ряд сходящийся = в долгосрочном периоде после появления шока происходит возврат к равновесию.


Остальные коэфы не интерпретируются, точнее интерпретируются с помощью графиков импульс респонс фанкщн 
Кроме альф значимы только краткосрочные шоки по переменной e для e = 0.7980*** и для U = -0.6200***.  

```{r}
plot(irf(VECmodel_const,impulse = "e",ci=0.95, boot =TRUE, ortho=FALSE))

plot(irf(VECmodel_const,impulse = "rw",ci=0.95, boot =TRUE, ortho=FALSE))

plot(irf(VECmodel_const,impulse = "U",ci=0.95, boot =TRUE, ortho=FALSE))
```
Долгосросчное равновесие может быть представлено в виде уравнения: e = -0.16 * rw +  3.94 * U. Параметр альфа для уравнения rw  = -0.1309*** (значим на 1% уровне) < 0, следовательно ряд сходящийся = в долгосрочном периоде после появления шока происходит возврат к равновесию. Параметр альфа для уравнения U  = -0.0312*** (значим на 1% уровне) < 0, следовательно ряд сходящийся = в долгосрочном периоде после появления шока происходит возврат к равновесию. Интерпретация краткосрочных шоков с помощью IRF: Если происходит краткосрочный шок по переменной e, то он будет влиять на саму переменную e положительно с самого начала и на протяжении минимум 10 лет, кроме того, на протяжении первых пяти лет со временем эффект нарастает, затем практически не меняется.  На переменную rw это оказывает отрицательное влияние в интервале от 1,5 до 3 лет с момента появления шока. На переменную U это оказывает отрицательное влияние с самого начала и на протяжении минимум 10 лет, кроме того, первые пять лет эффект усиливается, затем практически не меняется.  Если происходит краткосрочный шок по переменной rw, то как он будет влиять на переменную e на протяжении 10 лет сказать нельзя, так как ноль попадает в доверительный интервал, следовательно значения не значимы и не могут интерпретироваться. На саму переменную rw это оказывает стабильное положительное влияние с самого начала и на протяжении минимум 10 лет. На переменную U влияние так же не может быть проинтерпретировано так как ноль попадает в доверительный интервал, следовательно значения не значимы и не могут интерпретироваться.Если происходит краткосрочный шок по переменной U, то как он будет влиять на переменную e на протяжении 10 лет сказать нельзя, так как ноль попадает в доверительный интервал, следовательно значения не значимы и не могут интерпретироваться.  На переменную rw это оказывает отрицательное значимое влияние с нарастающим эффектом, начиная со второго года и минимум до конца 10, первые два года не могут быть проинтерпретированы. На саму переменную U это оказывает положительное значимое влияние со спадающим эффектов на протяжении первым двуз лет точно, затем влияние не может быть проинтерпретировано. 





```{r}
#VECM
library("tsDyn")
VECmodel_trend = VECM(M, lag = 1, r = 1, include = "trend",
estim = "ML")
summary(VECmodel_trend)
summary(VECmodel_trend)$coefMat
```

```{r}
plot(irf(VECmodel_trend,impulse = "e",ci=0.95, boot =TRUE, ortho=FALSE))

plot(irf(VECmodel_trend,impulse = "rw",ci=0.95, boot =TRUE, ortho=FALSE))

plot(irf(VECmodel_trend,impulse = "U",ci=0.95, boot =TRUE, ortho=FALSE))
```

```{r}
#VECM
library("tsDyn")
VECmodel_none = VECM(M, lag = 1, r = 1, include = "none",
estim = "ML")
summary(VECmodel_none)
summary(VECmodel_none)$coefMat
```

```{r}
plot(irf(VECmodel_none,impulse = "e",ci=0.95, boot =TRUE, ortho=FALSE))

plot(irf(VECmodel_none,impulse = "rw",ci=0.95, boot =TRUE, ortho=FALSE))

plot(irf(VECmodel_none,impulse = "U",ci=0.95, boot =TRUE, ortho=FALSE))
```

```{r}
#VECM
library("tsDyn")
VECmodel_both = VECM(M, lag = 1, r = 1, include = "both",
estim = "ML")
summary(VECmodel_both)
summary(VECmodel_both)$coefMat
```

```{r}
plot(irf(VECmodel_both,impulse = "e",ci=0.95, boot =TRUE, ortho=FALSE))

plot(irf(VECmodel_both,impulse = "rw",ci=0.95, boot =TRUE, ortho=FALSE))

plot(irf(VECmodel_both,impulse = "U",ci=0.95, boot =TRUE, ortho=FALSE))
```

```{r}
toLatex(summary(VECmodel), parenthese="Pvalue")
```

```{r}
BIC(VECmodel_const)
BIC(VECmodel_trend) #расходящийся
BIC(VECmodel_none) #расходящийся
BIC(VECmodel_both)
```
В ходе анализа были апробированы спецификации с включением константы (1) , тренда (2) , и константы, и тренда (3), и без включения константы и тренда (4). Вторая и четвертая спецификации были исключены, так как параметр альфа получался отрицательным и значимым, следовательно ряд являлся расходящимся = после возникновения шока не возвращается к равновесию в долгосрочной перспективе. Первая и третья спецификации сравнивались при помощи информационного критерия BIC. В первой спецификации соответствующее значение было меньше, что позволило сделать вывод о том, что спецификация с включением константы может быть названа лучшей для этих данных. Дополнительно для каждой модели строились ACF/PACF на остатках, по этому критерию все спецификации примерно равнозначны = не идеальны, но они приемлемы (совсем чуть-чуть вылезает).

```{r}
s1 = summary(VECmodel_const)
checkresiduals(s1$residuals[,1])
checkresiduals(s1$residuals[,2])
checkresiduals(s1$residuals[,3])
```


```{r}
s2 = summary(VECmodel_trend)
checkresiduals(s2$residuals[,1])
checkresiduals(s2$residuals[,2])
checkresiduals(s2$residuals[,3])
```


```{r}
s3 = summary(VECmodel_none)
checkresiduals(s3$residuals[,1])
checkresiduals(s3$residuals[,2])
checkresiduals(s3$residuals[,3])
```


```{r}
s4 = summary(VECmodel_both)
checkresiduals(s4$residuals[,1])
checkresiduals(s4$residuals[,2])
checkresiduals(s4$residuals[,3])
```


лучше еще тест на нормальность закинуть для каждой модели но в целом выбираем между моделями где альфа отрицательная и значимая, чтобы ряд сходился, и плюс по бику, затем у лучшей проверяем остатки на нормальность и с помощью чекрезидьюалс, вылезает ли что-то ( даже если вылезает, так честно и написать) затем интерпретируем импульс респонс фанкшн.

```{r}
plot(irf(VECmodel,impulse = "e",ci=0.95, boot =TRUE, ortho=FALSE))

plot(irf(VECmodel,impulse = "rw",ci=0.95, boot =TRUE, ortho=FALSE))

plot(irf(VECmodel,impulse = "U",ci=0.95, boot =TRUE, ortho=FALSE))
```