Include number of observations in textreg ARIMA output

Question

I'm running several time series regressions. My example is using different data, and I have taken NO care in making these models correct - they are just illustrative for my question. Sorry, I used a downloaded dataset rather than a created one. I hope that's okay.

My task is printing texreg outputs for various models, including but not limited to ARIMA forecasts. I need to include the number of observations for all models. The texreg extract functions do not include the number of observations for ARIMA objects - or at least the number of observations will never print in the gof statistics.

Two options I've thought of, and I'm hoping for help implementing either or both of them:

Add a custom row to the gof statistics with the number of observations included. This should be possible using custom.gof.rows, as discussed here and here. However, custom.gof.rows does not appear to a part of the texreg package as of version 1.36.23. Perhaps it was taken out for some reason? Not sure.

If there is a secret version of texreg that includes custom.gof.rows, can anyone link to it?

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))

This is the output, which I had to link (sorry)

Note that the model objects themselves aren't the reason (I think). I can extract the number of observations myself and print them separately.

#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))

Modify the extract function to include the number of observations so that they will automatically be included, as discussed the extended documentation here. This I'm completely new to, but I've made an attempt below to modify the textreg extract.Arima function to retrieve the number of observations using 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:

This runs but breaks the textreg function. If I'm on the wrong track, I'd also be happy to be linked to a tutorial for changing the extract functions generally. You may also notice that I'm using coeftest() in the code above to extract the Auto-ARIMA coefficients. I'm tangentially interested in writing an exract function for Auto-ARIMA forecast objects. I think this would be necessary anyway, because nobs only works on b itself, and not on coeftest(b). That's an aside though - I can figure that out when I get there.

Can anyone help with either or both avenues, or perhaps propose a different method for including the number of observations in a texreg output?

Thank you very much for any help.

Philip Leifeld Philip Leifeld · Accepted Answer · 2020-05-01T10:02:19

The custom.gof.rows argument was added to the texreg, screenreg etc. functions about a year ago and has not found its way to CRAN yet, but will eventually. For the time being, you can install the most recent version from GitHub if you'd like to use this argument. You can do so as follows (and will need to have installed the remotes package first):

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

The custom.gof.rows argument takes a named list of vectors. You were simply providing two values for the number of observations as separate arguments to the custom.gof.rows argument. Instead, you need to wrap both in a vector and attach a name to it, for example "Num. obs.". You can do this as follows:

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

That said, I have now added the number of observations to the extract methods for ARIMA models anyway, so you won't need this argument anymore.

Your attempt to modify the extract method on your own almost worked. You would have needed to include the include.nobs = TRUE argument also in the function head, and you would have needed to register the function as a method for the generic extract function as described in the article in the Journal of Statistical Software.

To make your life easier, I have done this for you. So if you install the most recent version of texreg from GitHub as per the instructions above, it should automatically include the number of observations.

The arima function you are using is defined in the stats package. It returns an Arima object. I have updated its extract method accordingly. The auto.arima function you are using above is a wrapper for that same function and returns an object with class forecast_ARIMA (and formerly, and secondarily, ARIMA). This is a bit confusing, but it's important to be aware of these distinctions. I have updated both methods for you with the number of observations, so your code

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

should now return the following output:

======================================
                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

Note, however, that you are using the coeftest function from the lmtest package for model b, which has a separate class and hence extract method and does not return the number of observations. To get the number of observations also for Model b, you could include it directly, without coeftest:

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

This returns the following output:

=======================================================
                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

If you do insist on using coeftest, you can save the output of the extract function in an intermediate object and manipulate this object. In the following code, I take the GOF block from the version without coeftest and inject it into the version with coeftest and then hand over the manipulated object instead of the original one to the screenreg function:

tr_b_coeftest <- extract(coeftest(b))
tr_b_plain <- extract(b)
[email protected] <- [email protected]
tr_b_coeftest@gof <- tr_b_plain@gof
[email protected] <- [email protected]
screenreg(list(tr_b_coeftest, a), single.row = TRUE)

This returns the following output:

=======================================================
                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

Include number of observations in textreg ARIMA output

1 Answers