rjdverse / rjdemetra

R interface to JDemetra+ v 2.x
https://rjdverse.github.io/rjdemetra
51 stars 16 forks source link

Joint F- Test p-value for User-defined calendar variables #56

Closed ssmite closed 5 years ago

ssmite commented 5 years ago

Hello! Is there a way how to get Joint F- Test p-value for User-defined calendar variables in RJDemetra?

joint_ftest

AQLT commented 5 years ago

The joint F-test is not exported but you can export all the informations needed to compute it: variance-covariance matrix, coefficients and degrees of freedom. Here one example to compute the F-test using car::linearHypothesis().

library(RJDemetra)
library(car)
#> Le chargement a nécessité le package : carData
rmod <- x13(ipi_c_eu[,"FR"],userdefined = "preprocessing.model.covar")
rmod$regarima
#> y = regression model + arima (0, 1, 1, 0, 1, 1)
#> Log-transformation: no
#> Coefficients:
#>           Estimate Std. Error
#> Theta(1)   -0.5270      0.048
#> BTheta(1)  -0.4865      0.051
#> 
#>               Estimate Std. Error
#> Monday       -0.133839      0.164
#> Tuesday      -0.002384      0.163
#> Wednesday     0.241712      0.163
#> Thursday     -0.531275      0.163
#> Friday        0.432474      0.164
#> Saturday      0.152956      0.163
#> Leap year    -0.045977      0.501
#> Easter [1]   -1.094082      0.335
#> LS (11-2008) -8.441602      1.307
#> LS (1-2009)  -7.274012      1.306
#> LS (5-2008)  -5.020079      1.257
#> 
#> 
#> Residual standard error: 1.665 on 323 degrees of freedom
#> Log likelihood = -624.7, aic =  1277 aicc =  1279, bic(corrected for length) = 1.252
coef.SA <- function(object, ...){
    object$regarima$regression.coefficients[,"Estimate"]
}
vcov.SA <- function(object, ...){
    result <- object$user_defined$preprocessing.model.covar
    rownames(result) <- colnames(result) <- rownames(object$regarima$regression.coefficients)
    result
}
df.residual.SA <- function(object, ...){
    object$regarima$loglik["neffectiveobs",] - object$regarima$loglik["np",]
}
linearHypothesis(rmod, c("Monday", "Tuesday", "Wednesday", "Thursday", "Friday", "Saturday")
                 , c(0,0,0,0,0,0), test = "F")
#> Linear hypothesis test
#> 
#> Hypothesis:
#> Monday = 0
#> Tuesday = 0
#> Wednesday = 0
#> Thursday = 0
#> Friday = 0
#> Saturday = 0
#> 
#> Model 1: restricted model
#> Model 2: rmod
#> 
#>   Res.Df Df      F   Pr(>F)   
#> 1    315                      
#> 2    309  6 2.9924 0.007418 **
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Created on 2019-08-14 by the reprex package (v0.3.0)

ssmite commented 5 years ago

I tried but I got an error : Error in vcov.default(model, complete = FALSE) : there is no vcov() method for models of class regarima, TRAMO_SEATS

Does it works only for x13 method?

AQLT commented 5 years ago

The example above works fine changing x13 by tramoseats. However it doesn't work with regARIMA/TRAMO you have to adapt the code to write methods for regarima objects. If you need help see what is code for regarima objects here: https://github.com/AQLT/rjdemetra/blob/master/R/regarima_generic.R

AQLT commented 4 years ago

With the last pull request you can directly use the linearHypothesis function:

library(RJDemetra)
library(car)
rmod <- x13(ipi_c_eu[,"FR"])
rmod2 <- tramoseats(ipi_c_eu[,"FR"])
linearHypothesis(rmod, c("Monday", "Tuesday", "Wednesday", "Thursday", "Friday", "Saturday")
                 , c(0,0,0,0,0,0), test = "F")
#> Linear hypothesis test
#> 
#> Hypothesis:
#> Monday = 0
#> Tuesday = 0
#> Wednesday = 0
#> Thursday = 0
#> Friday = 0
#> Saturday = 0
#> 
#> Model 1: restricted model
#> Model 2: rmod
#> 
#>   Res.Df Df      F   Pr(>F)   
#> 1    315                      
#> 2    309  6 2.9924 0.007418 **
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
linearHypothesis(rmod2, c("Monday", "Tuesday", "Wednesday", "Thursday", "Friday", "Saturday")
                 , c(0,0,0,0,0,0), test = "F")
#> Linear hypothesis test
#> 
#> Hypothesis:
#> Monday = 0
#> Tuesday = 0
#> Wednesday = 0
#> Thursday = 0
#> Friday = 0
#> Saturday = 0
#> 
#> Model 1: restricted model
#> Model 2: rmod2
#> 
#>   Res.Df Df      F   Pr(>F)   
#> 1    314                      
#> 2    308  6 2.8795 0.009588 **
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Created on 2019-10-13 by the reprex package (v0.3.0)