main.Rmd

---
title: "TS | Apple | Copart inc. | VAR"
author: "Krasnukhina"
date: "2022-12-12"
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)
library(lubridate)

data = read.csv("data.csv")
data = as.xts(data[, -c(1, 2)], order.by = as.Date(data$date))
colnames(data) = c('Apple', 'Copart_inc')
is.null(data)
```

Выборка включает данные о скоректированных на выплаченные дивиденды и проведенные сплиты ценах закрытия акций двух компаний из индекса NASDAQ 100: Apple и COPART Inc. Данные по котировкам акций взяты за каждый день работы биржи с 17 марта 1994 года по 27 марта 2018 года, соответственно, данные не представлены по датам, связанным с выходными и праздничными днями. Однако, поскольку биржа не работала в эти дни, то соответствующие наблюдения нельзя считать пропущенными значениями.

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

Выбранные для исследования компании Apple и COPART Inc. относятся к высокотехнологичным компаниям и включаются в единый индекс - NASDAQ 100, в чем проявляется их отраслевое сходство, однако специфика производимых ими товаров и услуг  имеет отличия. Кроме того, компания Apple по праву причисляется к числу IT-гигантов, обладая наибольшей капитализацией, в то время, как COPART Inc., в свою очередь, находится в числе компаний с наименьшей капитализацией среди включенных в индекс NASDAQ 100. В контексте выявленных различий актуальным становится вопрос, влияет ли изменение доходности акций компании с большим размером капитализации на доходность компании с меньшим размером капитализации и наоборот. 
Кроме того, актуальность исследования определяется также отсутствием единого мнения в научной среде о длительности сохранения эффекта при его наличии, а также о наличии или отсутствии авторегрессионного эффекта во взаимовлиянии.
Бенифициарами данного исследования, в первую очередь, являются инвесторы, так как оно вносит вклад в анализ факторов трансформации рыночной стоимости компании и определение природы колебаний котировок ценных бумаг, что может позволить снизить риски при формировании портфеля и увеличить потенциальную доходность.

В первом столбце данных представлены цены акций компании Apple, во втором столбце компании COPART Inc. Построим графики по каждой переменной:

```{r}
autoplot(data)
```

Можно увидеть, что в данных имеется тренд и дисперсия не является постоянной. Соответственно, можно предположить, что обе переменные нестационарны. Проверим, так ли это при помощи PP, KPSS и ADF тестов:

```{r}
# Apply Phillips-Perron test
pp.test(data[,1])
```
p-value >0.99, следовательно нулевая гипотеза о нестационарности ряда не отклоняется на любом уровне значимости. Ряд нестационарен. 

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

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

```{r}
# Apply Phillips-Perron test
pp.test(data[,2])
```
p-value >0.99, следовательно нулевая гипотеза о нестационарности ряда не отклоняется на любом уровне значимости. Ряд нестационарен. 

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

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

Поскольку оба ряда являются нестационарными, необходимо их преобразовать. Так как ранее на графиках мы увидели, что в данных есть тренд и дисперсия непостоянна, то для избавления от этих недостатков необходимо логарифмировать данные(это поможет снизить разброс дисперсии) и взять первые разности для избавления от тренда. Стоит заметить, что данные преобразования также преобразуют наши данные по ценам в доходности, что позволит легче интерпретировать полученные результаты анализа. Формула доходности:
$$r_t = \ln\left(\frac{P_t}{P_{t-1}}\right)$$
Преобразование данных:
```{r}
log = log(data)
autoplot(log)
log = as.xts(log)

# First difference
diff_data = diff.xts(data, 1, 1)
#remove first NA value
diff_data = diff_data[-1, ]
autoplot(diff_data)

# First difference
log_diff = diff.xts(log, 1, 1)
#remove first NA value
log_diff = log_diff[-1, ]
autoplot(log_diff)
```

При проведении преобразований можно увидеть, что применение логарифмирования и первых разностей по отдельности не делает ряды стационарными, что логично, поскольку данные преобразования борятся с разными проблемами. Поэтому далее будут использоваться разности логарифмированных значений.Как уже было сказано ранее, данные преобразования приводят нас к работе не с самими ценами, а с доходностями компаний. Использование доходностей обосновывается не только тем, что изначальные данные нестационарны и тем, что результаты будет удобнее интерпретировать, а также и тем, что ...

Также можно увидеть, что в данных присутствуют выбросы. Избавимся от них:

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

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

Построим графики по данным после преобразований и удаления выбросов:
```{r}
autoplot(log_diff)
```

Проверим полученные данные на стационарность:

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

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

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

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

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

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

```{r}
ggAcf(log_diff$Apple) + labs(title = "ACF PLot of Apple returns")
ggAcf(log_diff$Copart_inc) + labs(title = "ACF PLot of Copart inc. returns")
```

Apple: 6, 12, 24, 26, 30 лаги выходят за пределы доверительного интервала.
Copart inc.: 1, 7, 10, 13, 14, 21, 28 лаги выходят за пределы доверительного интервала.
Согласно анализу корролелограмм для каждой из исследуемых переменных большое количество лагов выходит за пределы доверительного интервала, что говорит о наличии автокорреляции в изучаемых временных рядах. Однако сезонности в данных не наблюдается. Кроме того, поскольку у нас 2 переменных в данных и целью исследования является не прогноз, а описание влияния переменных друг на друга, необходимо моделировать векторную авторегрессию. Также, поскольку данные нестационарны, а разности логарифмов стационарны, мы выбираем между моделью VAR и VECM в зависимости от наличия коинтеграции. Количество лагов для спецификации в этом случае выбирается с помощью информационных критериев. Реализуем выбор с помощью функции VARselect (Примечание: в случае построения модели VECM необходимо будет взять на один лаг меньше, так как специфика модели подразумевает расчет первых разностей в левой части уравнения модели):


Для выбора модели обратим внимание на данные и цель исследования: поскольку у нас 2 переменных в данных и целью исследования является не прогноз, а описание влияния переменных друг на друга, а также поскольку данные нестационарны, а разности логарифмов стационарны, мы выбираем между моделью VAR и VECM. Для выбора модели выберем количество лагов и проверим данные на наличие коинтеграции:
```{r}
VARselect(log_diff, lag.max = 20, type = "none", season = NULL, exogen = NULL)
```
По результатам выбора между спецификаций моделей с количеством лагов от 1 до 20 по критериям HQ и SC наилучшим количеством лагов будет 1 лаг. Проверим данные на коинтеграцию. Поскольку для модели VECM берется количество лагов на единицу меньшее, чем для модели VAR, но по нашим данным для VAR оптимальное количество лагов равно 1, то для проверки на наличие коинтеграции зададим количество лагов в VAR равное 2, поскольку 1 даст 0 как количество лагов в VECM:
```{r}
#cointegration test
library(urca)
jotest=ca.jo(log, type="eigen", K=2,
ecdet="none", spec="longrun")
summary(jotest)
```
```{r}
#cointegration test
library(urca)
jotest=ca.jo(log, type="trace", K=2,
ecdet="none", spec="longrun")
summary(jotest)
```
Для определения ранга коинтеграции были проведены виды тестов Йохансена: Trace и Eigen. При значении ранга 0 тестовая статистика равна 6.64 и 6.65 соответственно для Eigen и Trace тестов. Данные значения статистик меньше квантильных на уровнях значимости 10, 5 и 1 процент для обоих тестов для ранга равного 0. Нулевые гипотезы обоих тестов говорят о том, что значение ранга равно рассматриваемому r, альтернативная гипотеза для теста Eigen может быть сформулирована как rang = r + 1, для теста Trace альтернативная гипотеза rang > r. Нулевая гипотеза обоих тестов не отклоняется на 1% уровне значимости для r = 0, значит ранг коинтеграции равен 0 и коинтеграция отсутствует.

Мы проверяем наличие коинтеграции на данных после логарифмирования, так как в ходе построения модели VECM автоматически рассчитываются первые разности для левой части уравнения модели. Дополнительно мы проверяли наличие коинтеграции на исходных данных, была выявлена коинтеграция с рангом 1, однако если исходные данные подавать в модель, то после расчета первых разностей в левой части уравнения ряды не будут стационарными. Преобразования в виде первых разностей недостаточно для обеспечения стационарности, так как дисперсия данных  меняется со временем (график в начале файла, это проверялось с помощью (ADF, PP, KPSS тестов). Так как коинтеграция на прологарифмированных данных отсутствует, необходимо строить VAR in differences. В нашем случае это равноценно построению модели VAR на данных после перехода от котировок скорректированных цен закрытия к доходностям (log-diff преобразование) c одним лагом.

Таким образом, поскольку коинтеграция отсутствует, мы выбираем модель VAR с одним лагом и строим ее на логарифмических разностях изначальных данных. Система уравнений будет следующей:

$$Y_{1t} = {alpha_{1}} + {alpha_{11}}{Y_{1,t-1}} + {alpha_{12}}{Y_{2,t-1}} + {e_{1,t}}$$
$$Y_{2t} = {alpha_{2}} + {alpha_{21}}{Y_{2,t-1}} + {alpha_{22}}{Y_{2,t-1}} + {e_{2,t}}$$

```{r}
m2 = VAR(log_diff, p = 1, type = "const")
summary(m2)
```
По результатам построения модели мы получаем следующую систему уравнений:
$$A_{t} = {0.001} + {0.014}{C_{t-1}} + {0.002}{A_{t-1}} + {e_{1,t}}$$
$$C_{t} = {0.001} + {0.014}{A_{t-1}} - {0.047}{C_{t-1}} + {e_{2,t}}$$

Проверим чистоту остатков:

```{r}
ggAcf(residuals(m2))+ labs(title = "ACF PLot of model residuals")
```

По ACF графикам можно увидеть, что по обеим переменным некоторые лаги выходят за пределы доверительного интервала. Что говорит о том, что модель, к сожалению, учитывает не всю автокорреляцию, присутствующую в данных. Однако при выборе оптимального количества лагов у нас было два  варианта: 1 и 13, мы проверили дополнительно спецификацию с 13 лагами, она не позволяет учесть всю автокорреляцию, в графиках ACF по остаткам так же присутствуют лаги, выходящие за пределы доверительного интервала, однако значение информационного критерия BIC (исследование носит описательный, а не прогностический характер, поэтому используется именно этот критерий для принятия решения) увеличивается, что говорит о том, что модель с одним лагом и соответствующие остатки - не идеальный вариант, но наиболее оптимальных в ограничениях исследуемой выборки.

```{r}
roots(m2, modulus = TRUE)
```
Полученные значения < 1, что говорит о стационарности остатков, наша модель является приемлемой.

```{r}
checkresiduals(m2$varresult$Apple$residuals)
checkresiduals(m2$varresult$Copart_inc$residuals)
```

Кроме того, на графиках выше видно, что среднее примерно равно 0, что приближает наши остатки к виду white noise (требование модели VAR), однако заметно, что дисперсия не является константой, хотя и не изменяется драматично. О неравенстве дисперсии константе говорит также результат тестирования на гетероскедастичность остатков:

```{r}
arch.test(m2)
```
Очень маленькое значение p-value не превышает любого адекватного уровня значимости, из чего следует, что мы не может отвергнуть нулевую гипотезу теста о гомоскедастичности ряда, дисперсия остатков действительно не является константой.

```{r}
s1 = stability(m2, type= "OLS-CUSUM")
plot(s1)
```

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

```{r}
plot(irf(m2,ci=0.95, boot =TRUE, ortho=TRUE, n.ahead = 3, runs = 100))
```

По результатам построения Impulse Response Function, можно сделать следующие выводы:
1. Увеличение доходности акций Apple на 1 стандартное отклонение окажет значимое положительное влияние на доходность компании COPART Inc. продолжительностью в 1 период, далее 0 попадает в доверительный интервал, поэтому мы не можем оценить влияние.
2. Увеличение доходности акций COPART Inc. не окажет значимого влияния на доходность акций Apple.
3. Эффект от шока по обеим переменным на свои доходности акций будет сохранятся на протяжении одного периода до возвращения к равновесию.

Так, мы получили, что изменение доходности акций компании с большей капитализацией оказывает воздействие на доходность акций компании с меньшей капитализацией, но не наоборот. Доходность акций большой компании не зависит от доходности акций компании меньшего размера.

Для проверки результатов на устойчивость увеличим количество итераций до 1000:

```{r}
plot(irf(m2,ci=0.95, boot =TRUE, ortho=TRUE, n.ahead = 3, runs = 1000))
```

При увеличении количества итераций результаты остались практически такими же, таким образом, результаты являются устойчивыми.

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