Скачайте данные отсюда.

Данные

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)

Наблюдения и планы:

  1. У всех есть как минимум два режима.
  2. Канада ведет себя иначе, чем другие. Почему?
  3. В Австралии есть два различных режима (нормальные смеси). Почему?
  4. Мы можем оценить средние значения 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

Упражнение: примените его к другим странам, а также → Австралии, США, Канаде.

Всего наилучшего.