На этой неделе я буду прогнозировать рейтинги одобрения президента с помощью обобщенной аддитивной модели (GAM). Рейтинги одобрения президента следуют четкой сезонной схеме в течение восьми лет, в течение которых президент может занимать свой пост. В начале первого срока президенту предстоят общенациональные выборы, и рейтинги одобрения будут очень высокими. Однако, поскольку президент уделяет основное внимание управлению, а не общественному мнению, рейтинги одобрения начнут падать. Таким образом, рейтинг одобрения президента будет снижаться в конце первого года его пребывания в должности и в течение всего второго года его правления. На третий год пребывания у власти президенту предстоит переизбрание, и ему будет уделяться больше внимания общественному мнению. Таким образом, рейтинги одобрения должны расти по мере приближения выборов. Если президент будет избираться на второй срок, рейтинги одобрения будут следовать той же тенденции, что и в первый срок. Они будут разрушаться в течение первых нескольких лет пребывания у власти, поскольку президент выделяет больше ресурсов на реализацию своей повестки дня, а не на государственную политику. Тем не менее, рейтинг одобрения вырастет за последний год пребывания у власти, потому что президент сосредоточится на укреплении своего наследия.

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

Первоисточником данных был Проект американского президентства. Данные, использованные в этом анализе, представляют собой временной ряд среднемесячных рейтингов одобрения президента с июля 1941 г. по июнь 2021 г. В ряду было несколько пробелов из-за нечастых опросов в начале периода. Я собрал дополнительные данные от Ропера, чтобы заполнить эти пробелы. Оставшиеся пробелы я заполнил с помощью кусочно-линейной интерполяции. Интерполяция выполнялась только внутри администраций, а не между администрациями, чтобы гарантировать достоверность оценок. Всего в выборке 15 разных президентов и 952 месяца.

rm(list=ls())
library(tidyverse)
library(xts)
library(zoo)
library(forecast)
library(mgcv)
library(lubridate)
df <- readRDS('PresidentApprovalData.rds')

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

df <- as.xts(df[,c("Approval","PopularVote",'Years.In.Office','rgdp.pc','cpi','unemp','Elected')],as.Date(df$Date))
plot(df$Approval)

На приведенном ниже рисунке показано среднее одобрение президента за годы их пребывания в должности, а также рейтинги одобрения президентом Байдена, Трампа и Обамы. График лучше иллюстрирует сезонность в рейтингах одобрения. Кривая среднего рейтинга одобрения имеет форму буквы «W», предсказанную гипотезой сезонности.

Интересно отметить силу эффекта переизбрания. Рейтинги одобрения Обамы и Трампа выросли прямо перед их вторыми выборами. Таким образом, президенты видят рост рейтингов прямо перед своими вторыми выборами независимо от того, побеждает ли президент на выборах или от текущего состояния страны.

df <- readRDS('PresidentApprovalData.rds')
df$Years.In.Office <- round(df$Years.In.Office,4)
Biden <- df[df$President=="Joseph Biden",c('Years.In.Office','Approval')]
names(Biden) <- c('Years.In.Office','Biden')
Trump <- df[df$President=="Donald Trump",c('Years.In.Office','Approval')]
names(Trump) <- c('Years.In.Office','Trump')
Obama <- df[df$President=="Barack Obama",c('Years.In.Office','Approval')]
names(Obama) <- c('Years.In.Office','Obama')
names <- c('Years.In.Office','Obama')
Average.PA <- df %>%
  dplyr::select(Approval,Years.In.Office) %>%
  dplyr::filter(Years.In.Office<=8) %>%
  group_by(Years.In.Office) %>%
  dplyr::summarise(Approval.mean=mean(Approval),
            Approval.sd=sd(Approval)) %>%
  ungroup()
Average.PA <- left_join(Average.PA,Biden,by='Years.In.Office')
Average.PA <- left_join(Average.PA,Trump,by='Years.In.Office')
Average.PA <- left_join(Average.PA,Obama,by='Years.In.Office')
head(Average.PA)
Average.PA %>%
dplyr::select(Years.In.Office,'Mean'=Approval.mean,Biden,Trump,Obama) %>%
  pivot_longer(-Years.In.Office) %>%
  ggplot(aes(Years.In.Office,value,color=name)) + geom_point() + geom_line() + theme_classic() + labs(title = "Average Presidential Approval and Time in Office")

Модель

Я буду использовать полупараметрическую GAM, чтобы отразить влияние времени пребывания в должности и экономических условий на рейтинги одобрения президента. Затем я буду использовать модель ARIMA, чтобы зафиксировать оставшуюся структуру автокорреляции, присутствующую в данных. Я буду использовать пять предикторов для прогнозирования одобрения президента, предназначенных для учета влияния результатов работы президента, текущего состояния страны и эффекта времени пребывания в должности, описанного выше. PopularVote – это доля голосов, полученных президентом на последних выборах. Президенты, которые не были избраны, имеют значение PopularVote, равное нулю. Данные для этой переменной были получены из американского проекта президентства. Year.In.Office — количество лет, в течение которых президент находится у власти. Влияние этой переменной должно быть нелинейным и следовать модели сезонности в форме буквы «W», описанной в предыдущем разделе. Поддержка президента должна быть очень высокой сразу после его избрания, ослабевать в течение первого, второго и третьего года, а затем возрастать прямо перед переизбранием. В случае переизбрания рейтинг одобрения должен снижаться на протяжении большей части второго срока, но повышаться непосредственно перед уходом с поста. Elected – это фиктивная переменная, указывающая, был ли президент избран или назначен.

Наконец, я включил три экономические переменные; реальный ВВП на душу населения (rgdp.pc), инфляция (cpi) и безработица (unemp). Все три эти переменные отстают на один год и измеряются с шагом в один год. Данные по всем трем переменным были получены из FRED. Полный набор очищенных данных для этого анализа также можно найти на моем Github.

df <- as.xts(df[,c("Approval","PopularVote",'Years.In.Office','rgdp.pc','cpi','unemp','Elected')],as.Date(df$Date))
head(df)

GAM представляют собой тип полупараметрической модели статистического обучения. В GAM функциональная форма может быть указана напрямую путем включения линейных и полиномиальных членов или изучена из данных с использованием базисных преобразований. Базовые преобразования, которые я буду использовать, — это штрафованные кубические сплайны. Я буду моделировать PopularVote и Elected как линейную зависимость с одобрением президента. Years.In.Office, rgdp.pc, unemp и cpi будут преобразованы с использованием штрафных кубических сплайнов. Таким образом, модель выучит соответствующую функциональную форму, чтобы смоделировать влияние Years.In.Office на одобрение президента на основе данных.

Далее я зафиксирую автокорреляцию в рейтингах одобрения президента с помощью ARIMA. После оценки модели GAM я вычисляю остатки. Затем модель ARIMA подгоняется к остаткам GAM. Важно отметить, что стационарность здесь не имеет значения, потому что, установив сначала GAM, я сделал остатки стационарными (дополнительную информацию см. в следующем Посте Stack Exchange). Порядок модели ARIMA определяется путем нахождения модели, минимизирующей BIC, с использованием алгоритма пошагового поиска. Пошаговый поиск можно реализовать в R с помощью функции auto.arima. После подбора модели ARIMA я вычисляю подобранные значения ARIMA и повторно запускаю GAM с подобранными значениями ARIMA в качестве предиктора. Поскольку ARIMA оценивается по остаткам GAM, коэффициент ARIMA будет близок к единице.

Код, который я использую для оценки моей модели, показан ниже.

fit.gam.arima <- function(df) {
  GAM <- gam(Approval~s(Years.In.Office,bs = "cr",fx=T)+s(rgdp.pc,bs = "cr",fx=T)+s(unemp,bs = "cr",fx=T)+s(cpi,bs = "cr",fx=T)+
               PopularVote+Elected,data=df,family = gaussian)
  df$res <- residuals(GAM)
  arima <- auto.arima(df$res)
  df$arima.fit <- as.numeric(fitted(arima))
  GAM <- gam(Approval~s(Years.In.Office,bs = "cr",fx=T)+s(rgdp.pc,bs = "cr",fx=T)+s(unemp,bs = "cr",fx=T)+s(cpi,bs = "cr",fx=T)+
               PopularVote+Elected+arima.fit,data=df,family =gaussian)
  mod <- list(GAM,arima)
  names(mod) <- c('GAM','ARIMA')
  
  return(mod)
}

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

test <- df[round(0.75*nrow(df)):nrow(df),]
df <- df[1:round(0.75*nrow(df)),]

Результаты регрессии показаны ниже. ARIMA имеет три термина авторегрессии (AR) и один термин скользящего среднего (MA). В модели GAM статистически значимы все переменные, кроме выборов. Коэффициент ARIMA также очень близок к единице. Наконец, скорректированный квадрат R большой, со значением 0,886.

Mod <- fit.gam.arima(df)
summary(Mod$ARIMA)
summary(Mod$GAM)

На рисунке ниже показано влияние всех нелинейных переменных. Для преобразования пространства будет интерпретироваться только Years.In.Office. Прежде всего следует отметить, что ни одна из этих оценок не является причинно-следственной, и они просто определяют прогностические тенденции.

Имеются данные о W-образном эффекте сезонности. Результаты показывают, что президентская поддержка медленно снижается в течение первых трех лет их пребывания у власти. Затем на четвертом году происходит перегиб, когда действующему президенту грозит переизбрание, если президент отбывает второй срок, поддержка снижается до шестого года. Затем поддержка увеличивается в конце второго президентского срока. Обратите внимание, что в пример включен FDR, поэтому переменная Years.In.Office может быть больше восьми.

par(mfrow=c(2,2))
plot.gam(Mod$GAM)

На следующем рисунке показаны некоторые стандартные регрессионные методы диагностики. Гистограмма и график QQ остатков предполагает, что распределение остатков очень близко к стандартному нормальному, с немного более толстыми хвостами. Кроме того, точки лежат близко к линии под углом 45 градусов на графике Отклик против подогнанных значений. На последнем рисунке показаны временные ряды президентских рейтингов одобрения (черные) и значения, подобранные моделью (красные). Мы видим, что модель хорошо соответствует данным, поскольку подобранные значения близко совпадают с наблюдаемыми данными.

par(mfrow=c(2,2))
gam.check(Mod$GAM)
Plot.Trend <- function(df,mod) {
  plot(df$Approval)
  df$trend.line <- as.numeric(predict(mod))
  lines(df$trend.line,col='red')
  addLegend(legend.loc = "bottomleft", legend.names = c("Approval","Forecast"), 
            pch=15,col = c("black","red"))
}
par(mfrow=c(1,1))
Plot.Trend(df,Mod$GAM)

Моделирование движения вперед

Чтобы оценить предсказательную силу вне выборки, я буду использовать моделирование вперед. Моделирование вперед — это способ узнать, насколько хорошо модель справилась бы с прогнозированием рейтингов одобрения президента, если бы я использовал ее в течение последних 238 месяцев. В моделировании, начиная с первого месяца тестовой выборки (август 2001 г.), я буду использовать модель для прогнозирования рейтингов одобрения президента на следующие 12 месяцев. Затем я переделаю модель с первым наблюдением в наборе тестов, включенным в данные. Я буду использовать эту новую модель для прогнозирования на следующие 12 месяцев. Этот процесс продолжается для каждого наблюдения в наборе тестов.

predict.gam.arima <- function(df,mod,k=1,sd.cap=F,lb,ub,se=F) {
  arima.forecast <- as.numeric(forecast(mod$ARIMA,k)[["mean"]])
  df$arima.fit <- as.numeric(arima.forecast)
  preds <- predict(mod$GAM,df)
  if (se) {
    df$arima.fit <- as.numeric(forecast(mod$ARIMA, k)[["lower"]][,"95%"])
    lower <- predict(mod$GAM,df,se.fit=T)
    lower <- lower$fit - 1.96*lower$se.fit
    
    df$arima.fit <- as.numeric(forecast(mod$ARIMA, k)[["upper"]][,"95%"])
    upper <- predict(mod$GAM,df,se.fit=T)
    upper <- upper$fit+1.96*upper$se.fit
    
    preds <- list(preds,lower,upper)
    names(preds) <- c('fit','lower','upper')
    if (sd.cap) {
      preds$fit <- pmax(preds$fit,lb)
      preds$fit <- pmin(preds$fit,ub)
    }
  }
  if (sd.cap&!se) {
    preds <- pmax(preds,lb)
    preds <- pmin(preds,ub)
  }
  
  return(preds)
}
Pres.Approval.BT <- function(df,test,p=12,Constrain.Prediction,m=1.5) {
  Mod <- fit.gam.arima(df)
  dates <- index(test)
  wf.predictions <- NULL
  wf.df <- df[,c("Approval","Years.In.Office","rgdp.pc","cpi","unemp","Elected","PopularVote")]
for (i in 2:length(dates)) {
    if (length(dates)-i>=(p-1)) {
      n <- p-1
      pred.dat <- as.data.frame(test[i:c(i+n),c('PopularVote','Years.In.Office','rgdp.pc','cpi','unemp','Elected')])
    } else {
      n <- length(dates) - i
      pred.dat <- as.data.frame(test[i:c(i+n),c('PopularVote','Years.In.Office','rgdp.pc','cpi','unemp','Elected')])
    }
    pred.dat$rgdp.pc <- pred.dat$rgdp.pc[1]# on date i will not know what gdp was for date i + 12, only gdp for days 1:i
    pred.dat$cpi <- pred.dat$cpi[1]
    pred.dat$unemp <- pred.dat$unemp[1]
    pred.dat$PopularVote <- pred.dat$PopularVote[1] # will not know popular vote results before it happend 
    if (Constrain.Prediction) {
      ub <- as.numeric(wf.df$Approval[nrow(wf.df)])+m*sd(wf.df$Approval)
      lb <- as.numeric(wf.df$Approval[nrow(wf.df)])-m*sd(wf.df$Approval)
      preds <- predict.gam.arima(pred.dat,Mod,k=c(1+n),sd.cap = T,
                                 lb=max(lb,0),ub=max(ub,100))
    } else {
      preds <- predict.gam.arima(pred.dat,Mod,k=c(1+n),sd.cap = F)
    }
    if (length(preds)<p) {
      preds <- c(preds,rep(NA,p-c(n+1)))
    }
    names(preds) <- paste0('fc.',seq(1:p))
    preds <- as.xts(t(preds),dates[i])
    preds$wf.mean <- mean(wf.df$Approval)
    wf.predictions <- rbind(wf.predictions,preds)
    wf.df <- rbind(wf.df,test[i,c("Approval","Years.In.Office","rgdp.pc","cpi","unemp","Elected","PopularVote")])
    Mod <- fit.gam.arima(wf.df)
    if (i%%25==0|i==length(dates)) {
      cat(paste0(round(100*i/length(dates)),'% Complete. \n'))
    }
  }
  
  return(wf.predictions)
}
wf.predictions <- Pres.Approval.BT(df,test,p=12,Constrain.Prediction = F)

Чтобы получить доступ к производительности модели в моделировании шага вперед, мне нужен эталон для сравнения прогнозов. Ориентиром, который я использовал, является среднее значение рейтинга одобрения президента. Интересующая статистика представляет собой вне выборки R в квадрате. Квадрат R вне выборки сравнивает ошибку прогнозирования модели с ошибкой, если бы мы использовали среднее значение шага вперед для прогнозирования одобрения президента. Положительные значения указывают на то, что модель превосходит среднее значение, а отрицательные значения указывают на то, что модель уступает среднему значению. Я также рассмотрю среднеквадратичную ошибку (RMSE), чтобы оценить ошибку предсказания модели и графики предсказаний модели по сравнению с наблюдаемыми данными.

Pres.Aproval.BT.Res <- function(BT.Res,test) {
  OOSR2 <- NULL
  RMSE <- NULL
  for (i in 1:c(ncol(BT.Res)-1)) {
    pred <- na.omit(as.numeric(BT.Res[,paste0('fc.',i)]))
    pred <- as.xts(pred,index(test)[c(i+1):nrow(test)])
    tmp <- test$Approval
    tmp <- na.omit(cbind(tmp,pred,BT.Res$wf.mean))
    plot(tmp,main=paste(i,'Month Forecast'),col=c("black","red","green"))
    print(addLegend(legend.loc = "bottomleft", legend.names = c("Approval","Forecast","WF.Mean"), 
                    pch=15,col = c("black","red","green")))
    tmp2 <- sqrt(mean((tmp$Approval-tmp$pred)^2))
    tmp <- 1 - (sum((tmp$Approval-tmp$pred)^2)/sum((tmp$Approval-tmp$wf.mean)^2))
    OOSR2 <- c(OOSR2,tmp)
    RMSE <- c(RMSE,tmp2)
  }
  
  print(plot(1:c(ncol(BT.Res)-1),OOSR2))
  print(abline(0,0))
  print(plot(1:c(ncol(BT.Res)-1),RMSE))
  
  out <- data.frame(1:c(ncol(BT.Res)-1),OOSR2,RMSE)
  names(out) <- c('FC.Period','OOSR2','RMSE')
  return(out)
}
Pres.Aproval.BT.Res(wf.predictions,test)

Результаты демонстрируют, что модель обладает предсказательной силой; однако качество прогноза ухудшается по мере увеличения периода прогноза. При прогнозировании на один месяц модель имеет значение R вне выборки, равное 0,8. Однако квадрат R вне выборки уменьшается по мере увеличения длины прогноза и становится отрицательным при продолжительности прогноза 9 месяцев. Аналогичная тенденция наблюдается и с RMSE.

Просмотр графика прогноза, наложенного на временной ряд утверждения Президентом, также информативен. Графики подтверждают результаты R в квадрате вне выборки; точность прогноза снижается по мере увеличения длины прогноза. Однако они также демонстрируют, что модель иногда предсказывает большие скачки одобрения, когда их не было. Мы можем дополнительно повысить точность прогноза, наложив ограничения на ход прогноза. Я повторно запустил симуляцию вперед с ограничением, согласно которому прогноз рейтинга одобрения президента не может превышать плюс-минус 1,5 стандартных отклонения от последнего наблюдаемого рейтинга одобрения. Это ограничение на изменения прогноза должно предотвратить большие скачки, наблюдаемые в предыдущем моделировании.

wf.predictions <- Pres.Approval.BT(df,test,p=12,Constrain.Prediction=T,m=1.5)
Pres.Aproval.BT.Res(wf.predictions,test)

Наложение ограничений плюс-минус 1,5 стандартного отклонения немного улучшает прогноз на один месяц. Квадрат R вне выборки увеличивается с 0,84 до 0,85. Однако мы видим гораздо больший прирост производительности для прогнозов на более длительные периоды. Квадрат R вне выборки для семимесячного прогноза увеличился с 0,17 до 0,18, а квадрат R вне выборки для восьмимесячного прогноза вырос с 0,05 до 0,08.

Прогноз одобрения Байдена

Результаты, представленные в предыдущем разделе, предполагают, что мы можем использовать модель для прогнозирования рейтингов одобрения президента на срок до восьми месяцев. Через восемь месяцев прогноз не лучше, чем при использовании выборочного среднего рейтинга одобрения. Прогноз рейтингов одобрения Байдена на следующие шесть месяцев показан ниже. Обратите внимание, что мои данные заканчиваются в июне. Поэтому прогноз охватывает период с июля 2021 года по февраль 2020 года.

Прогнозы предполагают, что рейтинги одобрения президента Байдена будут постепенно снижаться в течение следующих восьми месяцев. Точечные оценки предсказывают, что рейтинги одобрения президента Байдена достигнут 50% в октябре 2021 года и упадут ниже 50% в декабре 2021 года. В течение следующих восьми месяцев я буду следить за ежемесячными рейтингами одобрения и сравнивать их с модельными прогнозами.

df <- readRDS('PresidentApprovalData.rds')
Biden <- df[df$President=="Joseph Biden",]
df <- as.xts(df[,c("Approval","PopularVote",'Years.In.Office','rgdp.pc','cpi','unemp','Elected')],as.Date(df$Date))
Years.In.Office <- seq(as.numeric(Biden$Years.In.Office[nrow(Biden)])+0.0834,by=0.0834,length.out = 8)
rgdp.pc <- rep(as.numeric(Biden$rgdp.pc[nrow(Biden)]),8)
cpi <- rep(as.numeric(Biden$cpi[nrow(Biden)]),8)
unemp <- rep(as.numeric(Biden$unemp[nrow(Biden)]),8)
Elected <- rep(as.numeric(Biden$Elected[nrow(Biden)]),8)
PopularVote <- rep(as.numeric(Biden$PopularVote[nrow(Biden)]),8)
pred <- data.frame(Years.In.Office,rgdp.pc,unemp,cpi,PopularVote,Elected)
Mod <- fit.gam.arima(df)
ub <- as.numeric(df$Approval[nrow(df)])+1.5*sd(df$Approval)
lb <- as.numeric(df$Approval[nrow(df)])-1.5*sd(df$Approval)
tmp <- predict.gam.arima(pred,Mod,k=8,sd.cap=T,lb=lb,ub=ub,se=T)
Biden <- Biden[,c('Date','Approval')]
Biden$lower <- Biden$Approval
Biden$upper <- Biden$Approval
Biden <- as.xts(Biden[,-1],Biden$Date)
pred$Approval <- tmp$fit
pred$lower <- tmp$lower
pred$upper <- tmp$upper
pred <- pred[,c('Approval','lower','upper')]
pred <- as.xts(pred,seq.Date(as.Date('2021-06-01'),length.out = nrow(pred)+1,by='month')[2:9])
Biden <- rbind(Biden,pred)
plot(Biden,col = c("black","red","red"),lty=c(1,2,2))
Biden