Загружаем библиотеку quantmod, которая позволяет быстро скачивать данные по цене акций, курсам валют и строить всякие графики для технического анализа
library(quantmod)
library(knitr)
opts_chunk$set(cache = TRUE, dev = "png", warning = FALSE)
# cache=TRUE --- включает кэширование для быстрого повторного запуска
# программы dev='png' переводит все графики в растровый формат png, там
# почти исчезают проблемы с кириллицей warning=FALSE убирает вывод на
# экран предупреждений
Загружаем данные по акциям Гугла:
getSymbols(Symbols = "GOOG", from = "2012-01-01", to = "2013-03-01")
## As of 0.4-0, 'getSymbols' uses env=parent.frame() and auto.assign=TRUE by
## default.
##
## This behavior will be phased out in 0.5-0 when the call will default to
## use auto.assign=FALSE. getOption("getSymbols.env") and
## getOptions("getSymbols.auto.assign") are now checked for alternate
## defaults
##
## This message is shown once per session and may be disabled by setting
## options("getSymbols.warning4.0"=FALSE). See ?getSymbol for more details
## [1] "GOOG"
Смотрим, что там загрузилось:
head(GOOG)
## GOOG.Open GOOG.High GOOG.Low GOOG.Close GOOG.Volume
## 2012-01-03 652.9 668.1 652.4 665.4 3676500
## 2012-01-04 665.0 670.2 660.6 668.3 2864000
## 2012-01-05 662.1 664.0 656.2 659.0 3282900
## 2012-01-06 659.1 660.0 649.8 650.0 2692900
## 2012-01-09 646.5 647.0 621.2 622.5 5822600
## 2012-01-10 629.8 633.8 616.9 623.1 4395600
## GOOG.Adjusted
## 2012-01-03 665.4
## 2012-01-04 668.3
## 2012-01-05 659.0
## 2012-01-06 650.0
## 2012-01-09 622.5
## 2012-01-10 623.1
Для удобства называем скорректированную на дивидент цену одной буквой
y <- GOOG$GOOG.Adjusted
Строим график
plot(y)
Даже по графику видно, что ряд y не стационарный. Например, если мысленно поделить график на две части, то среднее в левой явно не равно среднему в правой.
Такой график характерен для случайного блуждания, т.е. процесса примерно вида \( y_t=y_{t-1}+u_t \).
Для нестационарного ряда не имеет смысла расчёт ACF и PACF, так как корреляция между случайными величинами отстоящими на \( k \) шагов может меняться во времени. Тем не менее, глянем на графики:
Пара бессмысленных и беспощадных графиков
acf(y)
pacf(y)
Несмотря на свою бессмысленность, они полезны тем, что по ним можно опознать случайное блуждание. Если мы тупо попросили компьютер построить графики ACF и PACF и оказалось, что ACF крайне медленно убывает, а PACF обрывается в ноль после первой частной корреляции, то, вероятно, перед нами случайное блуждание.
Переходим к разностям, \( z_t=y_t - y_{t-1} \)
z = diff(y)
График для \( z \) уже похоже на стационарный, \( z \) крутитися вокруг среднего в полосе примерно постоянной ширины.
plot(z)
Графики ACF и PACF для \( z \)
acf(z, na.action = na.omit)
pacf(z, na.action = na.omit)
Оцениваем модель arima(1,1,0), arima(0,1,1), arima(1,1,1)
model1 <- arima(y, order = c(1, 1, 0))
model2 <- arima(y, order = c(0, 1, 1))
model3 <- arima(y, order = c(1, 1, 1))
Смотрим на оценки ARIMA(1,1,1)
print(model3)
##
## Call:
## arima(x = y, order = c(1, 1, 1))
##
## Coefficients:
## ar1 ma1
## -0.985 1.000
## s.e. 0.017 0.011
##
## sigma^2 estimated as 88.6: log likelihood = -1062, aic = 2130
Видим, что коэффициенты значимы. К сожалению, они не очень хорошо интерпретируются.
Строим прогноз по модели ARIMA(1,1,1) на 5 шагов вперед
predict(model3, n.ahead = 5)
## $pred
## Time Series:
## Start = 292
## End = 296
## Frequency = 1
## [1] 806.8 806.2 806.8 806.2 806.7
##
## $se
## Time Series:
## Start = 292
## End = 296
## Frequency = 1
## [1] 9.422 13.410 16.390 18.965 21.178
Компьютер оценивает все модели минимизируя или максимизируя какую-нибудь функцию. При применении МНК минимизируется сумма квадратов остатков, в методе максимального правдподобия — максимизируется функция правдоподобия. Для оптимизации придумано много алгоритмов. Но бывает, что они не срабатывают на какой-нибудь простой функции…
Найдите минимум данной функции устно, не используя компьютер: \[ f(x_1,x_2,x_3)=0.01\cdot (x_1-0.5)^2+|x_1^2-x_2|+|x_1^2-x_3| \]
А теперь попробуем с компьютером:
f <- function(x) {
return(0.01 * (x[1] - 0.5)^2 + abs(x[1]^2 - x[2]) + abs(x[1]^2 - x[3]))
}
nlm(f, c(0, 0, 0))
## $minimum
## [1] 0.002476
##
## $estimate
## [1] 2.490e-03 6.591e-06 6.591e-06
##
## $gradient
## [1] -0.017562 0.003906 0.003906
##
## $code
## [1] 3
##
## $iterations
## [1] 4
Это не проблема R, это проблема алгоритмов поиска экстремума — им тяжело работать с недифференциируемой функцией, у которой скачком меняется скорость роста.