Analytical vcov for dynamic alm(). Issue #52

config-i1 · Mar 9, 2021 · ecadd1c · ecadd1c
1 parent 553bba1
commit ecadd1c
Show file tree

Hide file tree

Showing 3 changed files with 55 additions and 38 deletions.
diff --git a/NEWS b/NEWS
@@ -1,9 +1,10 @@
-greybox v0.6.8 (Release data: 2021-03-03)
+greybox v0.6.8 (Release data: 2021-03-09)
 ==============
 
 Changes:
 * MAE, MSE etc now have parameter "holdout" instead of "actual". This is needed for consistency purposes with measures() function.
 * qqplot for alm() with Poisson and Negative Binomial distributions.
+* alm() with distribution=c("dpois","dnorm") and ar>0 now relies on analytical covariance matrix.
 
 Bugfixes:
 * A fix for an exotic case of zoo + data.table, resulting in destruction of response variable.

diff --git a/R/alm.R b/R/alm.R
@@ -654,14 +654,14 @@ alm <- function(formula, data, subset, na.action,
     #### Define the rest of parameters ####
     ellipsis <- list(...);
     # If arima was provided in the old style
-    if(orders[1]==0 && !is.null(ellipsis$ar)){
-        orders[1] <- ellipsis$ar;
+    if(orders[1]==0 && !is.null(ellipsis[["ar", exact=TRUE]])){
+        orders[1] <- ellipsis[["ar", exact=TRUE]]$ar;
     }
-    if(orders[2]==0 && !is.null(ellipsis$i)){
-        orders[2] <- ellipsis$i;
+    if(orders[2]==0 && !is.null(ellipsis[["i", exact=TRUE]])){
+        orders[2] <- ellipsis[["i", exact=TRUE]];
     }
-    if(orders[3]==0 && !is.null(ellipsis$ma)){
-        orders[3] <- ellipsis$ma;
+    if(orders[3]==0 && !is.null(ellipsis[["ma", exact=TRUE]])){
+        orders[3] <- ellipsis[["ma", exact=TRUE]];
     }
 
     # Parameters for distributions
@@ -971,27 +971,6 @@ alm <- function(formula, data, subset, na.action,
         recursiveModel <- TRUE;
     }
 
-    #### Define what to do with the maxeval ####
-    if(is.null(ellipsis$maxeval)){
-        if(any(distribution==c("dchisq","dpois","dnbinom","dbcnorm","plogis","pnorm")) || recursiveModel){
-            maxeval <- 500;
-        }
-        # The following ones don't really need the estimation. This is for consistency only
-        else if(any(distribution==c("dnorm","dlnorm","dlogitnorm")) & !recursiveModel && any(loss==c("likelihood","MSE"))){
-            maxeval <- 1;
-        }
-        else{
-            maxeval <- 200;
-        }
-        # LASSO / RIDGE need more iterations to converge
-        if(any(loss==c("LASSO","RIDGE"))){
-            maxeval <- 1000;
-        }
-    }
-    else{
-        maxeval <- ellipsis$maxeval;
-    }
-
     dataWork <- eval(mf, parent.frame());
     dataTerms <- terms(dataWork);
     # Make this numeric, to address potential issues with zoo + data.table
@@ -1505,6 +1484,30 @@ alm <- function(formula, data, subset, na.action,
             denominator <- NULL;
         }
 
+        #### Define what to do with the maxeval ####
+        if(is.null(ellipsis$maxeval)){
+            if(any(distribution==c("dchisq","dpois","dnbinom","dbcnorm","plogis","pnorm")) || recursiveModel){
+                # maxeval <- 500;
+                maxeval <- length(B) * 40;
+            }
+            # The following ones don't really need the estimation. This is for consistency only
+            else if(any(distribution==c("dnorm","dlnorm","dlogitnorm")) & !recursiveModel && any(loss==c("likelihood","MSE"))){
+                maxeval <- 1;
+            }
+            else{
+                # maxeval <- 200;
+                maxeval <- length(B) * 40;
+            }
+            # LASSO / RIDGE need more iterations to converge
+            if(any(loss==c("LASSO","RIDGE"))){
+                # maxeval <- 1000;
+                maxeval <- length(B) * 80;
+            }
+        }
+        else{
+            maxeval <- ellipsis$maxeval;
+        }
+
         # Although this is not needed in case of distribution="dnorm", we do that in a way, for the code consistency purposes
         res <- nloptr(B, CF,
                       opts=list(algorithm=algorithm, xtol_rel=xtol_rel, maxeval=maxeval, print_level=print_level,

diff --git a/R/methods.R b/R/methods.R
@@ -748,15 +748,26 @@ vcov.alm <- function(object, bootstrap=FALSE, ...){
                            "It is recommended to use bootstrap=TRUE option in this case."),
                     call.=FALSE);
         }
-        arimaModel <- !is.null(object$other$arima);
 
-        if(!arimaModel && any(object$distribution==c("dnorm","dlnorm","dbcnorm","dlogitnorm"))){
-            matrixXreg <- object$data[,-1,drop=FALSE];
+        # If there are ARIMA orders, define them.
+        if(!is.null(object$other$arima)){
+            arOrders <- object$other$orders[1];
+            iOrders <- object$other$orders[2];
+            maOrders <- object$other$orders[3];
+        }
+        else{
+            arOrders <- iOrders <- maOrders <- 0;
+        }
+
+        if(iOrders==0 && maOrders==0 && any(object$distribution==c("dnorm","dlnorm","dbcnorm","dlogitnorm"))){
+            matrixXreg <- object$data;
             if(interceptIsNeeded){
-                matrixXreg <- cbind(1,matrixXreg);
+                matrixXreg[,1] <- 1;
                 colnames(matrixXreg)[1] <- "(Intercept)";
             }
-            # colnames(matrixXreg) <- names(coef(object));
+            else{
+                matrixXreg <- matrixXreg[,-1,drop=FALSE];
+            }
             matrixXreg <- crossprod(matrixXreg);
             vcovMatrixTry <- try(chol2inv(chol(matrixXreg)), silent=TRUE);
             if(any(class(vcovMatrixTry)=="try-error")){
@@ -784,14 +795,16 @@ vcov.alm <- function(object, bootstrap=FALSE, ...){
             }
             rownames(vcov) <- colnames(vcov) <- variablesNames;
         }
-        else if(!arimaModel && object$distribution=="dpois"){
-            matrixXreg <- object$data[,-1,drop=FALSE];
-            obsInsample <- nobs(object);
+        else if(iOrders==0 && maOrders==0 && object$distribution=="dpois"){
+            matrixXreg <- object$data;
             if(interceptIsNeeded){
-                matrixXreg <- cbind(1,matrixXreg);
+                matrixXreg[,1] <- 1;
                 colnames(matrixXreg)[1] <- "(Intercept)";
             }
-            # colnames(matrixXreg) <- names(coef(object));
+            else{
+                matrixXreg <- matrixXreg[,-1,drop=FALSE];
+            }
+            obsInsample <- nobs(object);
             FIMatrix <- matrixXreg[1,] %*% t(matrixXreg[1,]) * object$mu[1];
             for(j in 2:obsInsample){
                 FIMatrix[] <- FIMatrix + matrixXreg[j,] %*% t(matrixXreg[j,]) * object$mu[j];