Скачайте данные отсюда.
Данные
head(data) dim(data) > 7247,11 #not good head visualization
Типы переменных
str(data)
Важные переменные (здравый смысл)
Непрерывные переменные
- PercentofBaseline
Категориальные переменные
- Дата
- Город (внутри штата)
- Государство (внутри страны)
- Страна (мы сосредоточим внимание на этом важном категориальном факторе)
- Название аэропорта (внутри города)
Наблюдения
Это данные временного ряда.
Мы не знаем, есть ли у нас информация обо всех аэропортах, городах четырех стран. Таким образом, движение PecentofBaseline имеет значение только во времени в каждом из указанных аэропортов, городов и стран.
Исследовательский анализ данных
Одномерный
Плотность ядра PercentofBaseline
ggplot(data, aes(x = PercentOfBaseline, fill = "red"))+geom_density(alpha = 0.2)
Двумерный
Непрерывный против категориального
Плотность ядра в процентах от базовой линии по стране
ggplot(data, aes(x = PercentOfBaseline, fill = Country))+geom_density(alpha = 0.2)
Наблюдения и планы:
- У всех есть как минимум два режима.
- Канада ведет себя иначе, чем другие. Почему?
- В Австралии есть два различных режима (нормальные смеси). Почему?
- Мы можем оценить средние значения MLE для этой нормальной смеси и проверить, подходит ли она для нормальной смеси.
Упражнение: вы можете сделать это применительно к каждому из аэропортов, городов и штатов.
PercentofBaseline против времени (по стране)
ggplot(data, aes(x = Date, y = PercentOfBaseline, group = Country$, colour = Country)) + geom_line() +facet_wrap(~ Country)
Понимание
- В Австралии процент от базового уровня со временем снизился (из-за более эффективной государственной политики).
- В Канаде всегда высокий процент от базового уровня, что приводит к разному поведению.
- Чили имеет растущий процент от базового уровня.
- В целом по Соединенным Штатам процент от исходного уровня колеблется и высок.
Упражнение: вы можете сделать это применительно к каждому из аэропортов, городов и штатов.
Моделирование
Чили
lyr) data_chile = data %>% filter(Country == "Chile") chile = data_chile[,c(2,5)] head(chile)
Преобразование в структуру данных временного ряда *
- 1,2,3,4,… как временной индекс без значений NA.
library(zoo) z <- read.zoo(chile, format = "%Y-%m-%d")#zoo series for dates time(z) <- seq_along(time(z))#sequential data as 1,2,3,4,... ts_chile = as.ts(z) #conversion into time series data structure head(ts_chile)
- Дата как временной индекс со значениями NA, которые необходимо удалить.
library(zoo) library(tseries)#for removing na from ts() z <- read.zoo(chile, format = "%Y-%m-%d")#zoo series for dates ts_chile = as.ts(z) ts_chile = na.remove(ts_chile)#remove na from time series head(ts_chile)
Если вы конвертируете здесь, это создаст странные числа. Я предпочитаю делать это как 1,2,3,4,… временного индекса.
Постройте временной ряд
plot(ts_chile)
Мы обсудим две части временного ряда:
- Разложение (тренд + сезон + ошибка)
- Прогнозирование
Разложение
decompose(ts_chile)
Прогноз
Модель (модель ARIMA)
library(forecast) model = auto.arima(ts_chile) model
Участок АКФ
acf(model$residuals, main = "Correlogram")
Участок PACF
pacf(model$residuals, main = "Partial Correlogram")
Ljung Box Test
Тест определяет, являются ли ошибки iid (то есть белым шумом) или за ними стоит что-то еще; являются ли автокорреляции для ошибок или остатков ненулевыми.
Box.test(model$residuals, type = "Ljung-Box")
Мы не отвергаем гипотезу о том, что модель хорошо подходит.
Нормальность остатков
hist(model$residuals, freq = FALSE) lines(density(model$residuals))
Прогноз
forecast = forecast(model,4)#4 = number of units you want to library(ggplot2) autoplot(forecast) #plot of the model
accuracy(forecast) #performance of the model
forecast #forecast values
Упражнение: примените его к другим странам, а также → Австралии, США, Канаде.
Всего наилучшего.