Включить количество наблюдений в текстовый вывод ARIMA

Question

Я управляю несколькими регрессиями временных рядов. Мой пример использует разные данные, и я не позаботился о том, чтобы эти модели были правильными - они просто иллюстрируют мой вопрос. Извините, я использовал загруженный набор данных, а не созданный. Я надеюсь, что все в порядке.

Моя задача - распечатать выходные данные texreg для различных моделей, включая, помимо прочего, прогнозы ARIMA. Мне нужно включить количество наблюдений для всех моделей. Функции извлечения texreg не включают количество наблюдений для объектов ARIMA - или, по крайней мере, количество наблюдений никогда не будет напечатано в статистике gof.

Два варианта, о которых я подумал, и я надеюсь на помощь в реализации одного или обоих из них:

1. Добавление пользовательской строки в статистику gof с количеством наблюдений включен. Это должно быть возможным с использованием custom.gof.rows, как обсуждено здесь и здесь . Однако custom.gof.rows не отображается в составе пакета texreg начиная с версии 1.36.23. Возможно, это было удалено по какой-то причине? Не уверен.

Если есть секретная версия texreg, которая включает custom.gof.rows, кто-нибудь может дать ссылку на нее?

library('ggplot2')
library('forecast')
library('tseries')
library('texreg')
library('lmtest')

#data can be downloaded from: https://archive.ics.uci.edu/ml/datasets/Bike+Sharing+Dataset

#Change to data path
path<-c("")

#Import data
daily_data = read.csv(paste(path,'day.csv',sep=""), header=TRUE, stringsAsFactors=FALSE)

#Mark a shift date
daily_data$shift<- ifelse(daily_data$instant>12,1,0)

#Create time series data
dt <- ts(daily_data$temp, frequency=12, start= c(2011, 1))

#Define input vars for ARIMA
(xvars <- as.matrix(daily_data[, c("instant", "shift")]))

#Basic ARIMA
a<- arima(dt, xreg=xvars)

#Auto-ARIMA
b<- auto.arima(dt, xreg=xvars)

##Where I want to include number of observations either automatically or using custom.gof.rows
screenreg(list(coeftest(b),a))

Это вывод, который мне пришлось связать (извините)

Обратите внимание, что сами объекты модели не являются причина (я думаю). Я могу самостоятельно извлечь количество наблюдений и распечатать их отдельно.

#Both models have number of observations in their objects that I can extract
nobs_a<-nobs(a)
nobs_b<-nobs(b)

cat("Number of obs ARIMA: ",nobs_a,"\n")

cat("Number of obs Auto-ARIMA: ",nobs_b,"\n")

##Custom.gof.rows doesn't appear to be in texreg version 1.36.23
screenreg((list(coeftest(b),a)), custom.gof.rows = list(nobs_a,nobs_b))

Измените функцию извлечения, включив в нее количество наблюдений, чтобы они автоматически включались, как обсуждалось в расширенной документации здесь . Я совершенно новичок в этом, но ниже я попытался изменить функцию textreg extract.Arima, чтобы получить количество наблюдений, используя nobs.

function (model, include.pvalues = FALSE, include.aic = TRUE, 
          include.loglik = TRUE, ...)  {   mask <- model$mask   nam <- names(model$coef)   
   co <- model$coef   sdev <-
    sqrt(diag(model$var.coef))   if (include.pvalues == TRUE) {
    t.rat <- rep(NA, length(mask))
    t.rat[mask] <- co[mask]/sdev
    pt <- 2 * pnorm(-abs(t.rat))
    setmp <- rep(NA, length(mask))
    setmp[mask] <- sdev   }   else {
    pt <- numeric()
    setmp <- sdev   }   gof <- numeric()   gof.names <- character()   gof.decimal <- logical()   if 
     (include.aic == TRUE) {
    aic <- AIC(model)
    gof <- c(gof, aic)
    gof.names <- c(gof.names, "AIC")
    gof.decimal <- c(gof.decimal, TRUE)   }   

    ##This is the section I added - ripped more or less intact from the extract.lm function         
    if (include.nobs == TRUE) {
    gof <- c(gof, nobs)
    gof.names <- c(gof.names, "Num. obs.")
    gof.decimal <- c(gof.decimal, FALSE)   }   

    if (include.loglik == TRUE) {
    lik <- model$loglik
    gof <- c(gof, lik)
    gof.names <- c(gof.names, "Log Likelihood")
    gof.decimal <- c(gof.decimal, TRUE)   }   
    tr <- createTexreg(coef.names = nam, coef = co, se = setmp, 
                     pvalues = pt, gof.names = gof.names, gof = gof, gof.decimal = gof.decimal)   > 
      return(tr) } <bytecode:

Это работает, но нарушает функцию textreg. Если я нахожусь на неправильном пути, я также был бы рад получить ссылку на учебник по изменению функций извлечения вообще. Вы также можете заметить, что я использую coeftest () в приведенном выше коде для извлечения коэффициентов Auto-ARIMA. Я заинтересован в написании функции извлечения для объектов прогноза Auto-ARIMA. Я думаю, что это было бы необходимо в любом случае, потому что nobs работает только на b, а не на coeftest (b). Хотя это в стороне - я могу понять это, когда доберусь туда.

Может кто-нибудь помочь с одним или обоими направлениями, или, возможно, предложить другой метод для включения количества наблюдений в вывод texreg?

Большое спасибо за любую помощь.

Answer 1

Аргумент custom.gof.rows был добавлен к texreg, screenreg et c. функционирует около года go и еще не нашел свой путь в CRAN, но в конечном итоге будет. В настоящее время вы можете установить самую последнюю версию из GitHub, если хотите использовать этот аргумент. Вы можете сделать это следующим образом (и сначала нужно будет установить пакет remotes):

library(remotes)
install_github("leifeld/texreg")

Аргумент custom.gof.rows принимает именованный список векторов. Вы просто указали два значения для числа наблюдений в качестве отдельных аргументов для аргумента custom.gof.rows. Вместо этого вам нужно обернуть оба вектора и прикрепить к нему имя, например «Num. Obs.». Вы можете сделать это следующим образом:

screenreg(list(coeftest(b), a),
          custom.gof.rows = list("Num. obs." = c(nobs(b), nobs(a))))

Тем не менее, теперь я добавил количество наблюдений к extract методам для моделей ARIMA, так что вам больше не нужен этот аргумент.

Ваша попытка изменить метод extract самостоятельно почти сработала. Вы должны были бы включить аргумент include.nobs = TRUE также в заголовок функции, и вам нужно было бы зарегистрировать функцию как метод для функции generi c extract, как описано в статье в журнале. Статистического программного обеспечения .

Чтобы сделать вашу жизнь проще, я сделал это для вас. Поэтому, если вы устанавливаете самую последнюю версию texreg из GitHub в соответствии с инструкциями выше, она должна автоматически включать количество наблюдений.

Используемая вами функция arima определена в stats пакет. Возвращает объект Arima. Я обновил его extract метод соответственно. Функция auto.arima, которую вы используете выше, является оболочкой для той же самой функции и возвращает объект с классом forecast_ARIMA (и ранее, и во вторую очередь, ARIMA). Это немного сбивает с толку, но важно знать об этих различиях. Я обновил для вас оба метода, указав количество наблюдений, поэтому ваш код

screenreg(list(coeftest(b), a))

теперь должен возвращать следующий вывод:

======================================
                Model 1    Model 2    
--------------------------------------
ar1              0.96 ***             
                (0.04)                
ar2             -0.30 ***             
                (0.05)                
ar3              0.13 *               
                (0.05)                
ar4              0.02                 
                (0.05)                
ar5              0.17 ***             
                (0.04)                
intercept        0.43 **      0.21 ***
                (0.13)       (0.05)   
instant          0.00         0.00 *  
                (0.00)       (0.00)   
shift            0.02         0.25 ***
                (0.05)       (0.05)   
--------------------------------------
AIC                        -440.26    
BIC                        -421.88    
Log Likelihood              224.13    
Num. obs.                   731       
======================================
*** p < 0.001; ** p < 0.01; * p < 0.05

Обратите внимание, что вы используете Функция coeftest из пакета lmtest для модели b, которая имеет отдельный класс и, следовательно, метод extract и не возвращает количество наблюдений. Чтобы получить число наблюдений также для модели b, вы можете включить ее напрямую, без coeftest:

screenreg(list(b, a), single.row = TRUE)

. Это вернет следующий вывод:

=======================================================
                Model 1              Model 2           
-------------------------------------------------------
ar1                 0.96 (0.04) ***                    
ar2                -0.30 (0.05) ***                    
ar3                 0.13 (0.05) *                      
ar4                 0.02 (0.05)                        
ar5                 0.17 (0.04) ***                    
intercept           0.43 (0.13) **      0.21 (0.05) ***
instant             0.00 (0.00)         0.00 (0.00) *  
shift               0.02 (0.05)         0.25 (0.05) ***
-------------------------------------------------------
AIC             -2183.71             -440.26           
AICc            -2183.46                               
BIC             -2142.36             -421.88           
Log Likelihood   1100.86              224.13           
Num. obs.         731                 731              
=======================================================
*** p < 0.001; ** p < 0.01; * p < 0.05

Если вы настаивайте на использовании coeftest, вы можете сохранить выходные данные функции extract в промежуточном объекте и манипулировать этим объектом. В следующем коде я беру блок GOF из версии без coeftest и внедряю его в версию с coeftest, а затем передаю управляемый объект вместо исходного объекту screenreg:

tr_b_coeftest <- extract(coeftest(b))
tr_b_plain <- extract(b)
tr_b_coeftest@gof.names <- tr_b_plain@gof.names
tr_b_coeftest@gof <- tr_b_plain@gof
tr_b_coeftest@gof.decimal <- tr_b_plain@gof.decimal
screenreg(list(tr_b_coeftest, a), single.row = TRUE)

Возвращает следующий вывод:

=======================================================
                Model 1              Model 2           
-------------------------------------------------------
ar1                 0.96 (0.04) ***                    
ar2                -0.30 (0.05) ***                    
ar3                 0.13 (0.05) *                      
ar4                 0.02 (0.05)                        
ar5                 0.17 (0.04) ***                    
intercept           0.43 (0.13) **      0.21 (0.05) ***
instant             0.00 (0.00)         0.00 (0.00) *  
shift               0.02 (0.05)         0.25 (0.05) ***
-------------------------------------------------------
AIC             -2183.71             -440.26           
AICc            -2183.46                               
BIC             -2142.36             -421.88           
Log Likelihood   1100.86              224.13           
Num. obs.         731                 731              
=======================================================
*** p < 0.001; ** p < 0.01; * p < 0.05