Skip to contents

x13 with PICKMDL and partial concurrent possibilities

x13_pickmdl.Rd

x13 can be run as usual (automdl) or with a PICKMDL specification. The ARIMA model, outliers and filters can be identified at a certain date and then held fixed (with a new outlier-span).

Usage

x13_pickmdl(
  series,
  spec,
  corona = FALSE,
  ...,
  pickmdl_method = "first",
  star = 1,
  when_star = warning,
  when_automdl = message,
  when_finalnotok = NULL,
  identification_end = NULL,
  identification_estimate.to = NULL,
  identify_t_filter = FALSE,
  identify_s_filter = FALSE,
  identify_outliers = TRUE,
  identify_arima_mu = TRUE,
  automdl.enabled = FALSE,
  fastfirst = TRUE,
  verbose = FALSE,
  output = "sa",
  add_comment = TRUE,
  old_crit2 = NA
)

x13_automdl(..., automdl.enabled = TRUE)

Arguments

series

x13 parameter

spec

An x13_spec output object or a list of several objects as outputted from x13_spec_pickmdl. In the case of a single object and when automdl.enabled is FALSE, spec will be converted internally by x13_spec_pickmdl with default five arima model specifications.

corona

Whether to update spec by outliers according to corona_outliers. FALSE or NULL means no update. TRUE or "ssb" means update.

...

Further x13 parameters (currently only parameter userdefined is additional parameter to x13).

pickmdl_method

crit_selection parameter or one of the two extra possibilities, "first_automdl" or "first_tryautomdl". In both cases the crit_selection parameter is "first" and the automdl model is added as the last pickmdl model.

  • "first_automdl": The automdl model is chosen whenever no pickmdl model is ok. In other words, the star parameter changes.

  • "first_tryautomdl": When no pickmdl model is ok: The automdl model is chosen if this model is ok, otherwise the star model is chosen.

star

crit_selection parameter

when_star

crit_selection parameter

when_automdl

Function to be called when automdl since no pickmdl model ok. Supply NULL to do nothing.

when_finalnotok

Function to be called, e.g. warning, when final run with final model is not ok. Supply NULL to do nothing. See crit_ok.

identification_end

To shorten the series before runs used to identify (arima) parameters. That is, the series is shortened by window(series, end = identification_end).

identification_estimate.to

To set x13_spec parameter estimate.to before runs used to identify (arima) parameters. This is an alternative to identification_end.

identify_t_filter

When TRUE, Henderson trend filter is identified by the shortened (see above) series.

identify_s_filter

When TRUE, Seasonal moving average filter is identified by the shortened series.

identify_outliers

When TRUE, Outliers are identified by the shortened series.

identify_arima_mu

When TRUE, arima.mu is identified by the shortened series (see arima_mu).

automdl.enabled

Logical value or any other value.

  • When set to FALSE (the default), the pickmdl routine is applied.

  • When set to TRUE, the automdl routine is performed.

  • For any value other than TRUE or FALSE, the ARIMA model is chosen as specified by spec.

Note that when automdl.enabled is not FALSE, if spec is a list containing several objects outputted from x13_spec_pickmdl, only the first object is used.

fastfirst

When TRUE and when pickmdl with crit_selection parameter "first", only as many models as needed are run. This affects the output when output = "all".

verbose

Printing information to console when TRUE.

output

One of "sa" (default), "spec" (final spec), "sa_spec" (both) and "all". See examples.

add_comment

When TRUE, a comment attribute (character vector with ok, ok_final and mdl_nr) will be added to the x13 output object. Use comment to get the attribute or ok to get the attribute converted to a list.

old_crit2

Logical. The p-value criterion used for PICKMDL criterion number 2. Set to FALSE for "Ljung-Box" and to TRUE for "Ljung-Box (residuals at seasonal lags)". This parameter can be overridden by setting the "pickmdl.old_crit2" option, in which case the option value will take precedence. The default value (NA) means that old_crit2 = (date_found < "2024-10-15"), where date_found refers to the end date according to identification_end or identification_estimate.to. If none of these is specified, old_crit2 is set to FALSE. The default value is chosen to ensure that the new criterion is phased in automatically.

Value

By default an x13 output object, or otherwise a list as specified by parameter output.

Examples

myseries <- pickmdl_data("myseries")

spec_a  <- x13_spec(spec = "RSA3", transform.function = "Log")

a <- x13_pickmdl(myseries, spec_a, verbose = TRUE)
#> old_crit2 set to FALSE since no identification specification. 
#>  p  d  q bp bd bq 
#>  2  1  0  0  1  1 
#> outlier.from updated: 2021-10-01  New outliers: 2011-06 
comment(a)
#>       ok ok_final   mdl_nr 
#>   "TRUE"   "TRUE"      "3" 
ok(a)
#> $ok
#> [1] TRUE
#> 
#> $ok_final
#> [1] TRUE
#> 
#> $mdl_nr
#> [1] 3
#> 
unlist(ok(a))
#>       ok ok_final   mdl_nr 
#>        1        1        3 
a$regarima
#> y = regression model + arima (2, 1, 0, 0, 1, 1)
#> Log-transformation: yes
#> Coefficients:
#>           Estimate Std. Error
#> Phi(1)      0.9848      0.061
#> Phi(2)      0.5766      0.061
#> BTheta(1)  -0.8921      0.045
#> 
#>             Estimate Std. Error
#> TC (6-2011)  -0.1905      0.047
#> 
#> 
#> Residual standard error: 0.06728 on 183 degrees of freedom
#> Log likelihood = 230.6, aic =  1260 aicc =  1260, bic(corrected for length) = -5.286
#> 

a2 <- x13_pickmdl(myseries, spec_a, identification_end = c(2014, 2))
a2$regarima
#> y = regression model + arima (0, 1, 1, 0, 1, 1)
#> Log-transformation: yes
#> Coefficients:
#>           Estimate Std. Error
#> Theta(1)   -0.8434      0.040
#> BTheta(1)  -0.9735      0.044
#> 
#> 
#> Residual standard error: 0.0707 on 185 degrees of freedom
#> Log likelihood =   216, aic =  1285 aicc =  1285, bic(corrected for length) = -5.243
#> 

# As above, another way
a3 <- x13_pickmdl(myseries, spec_a, identification_estimate.to = "2014-03-01")
a3$regarima
#> y = regression model + arima (0, 1, 1, 0, 1, 1)
#> Log-transformation: yes
#> Coefficients:
#>           Estimate Std. Error
#> Theta(1)   -0.8434      0.040
#> BTheta(1)  -0.9735      0.044
#> 
#> 
#> Residual standard error: 0.0707 on 185 degrees of freedom
#> Log likelihood =   216, aic =  1285 aicc =  1285, bic(corrected for length) = -5.243
#> 

a4 <- x13_automdl(myseries, spec_a, identification_end = c(2014, 2))
a4$regarima
#> y = regression model + arima (1, 0, 2, 0, 1, 1)
#> Log-transformation: yes
#> Coefficients:
#>           Estimate Std. Error
#> Phi(1)     -0.9522      0.030
#> Theta(1)   -1.1092      0.072
#> Theta(2)    0.3625      0.071
#> BTheta(1)  -1.0000      0.045
#> 
#>      Estimate Std. Error
#> Mean  0.02675      0.004
#> 
#> 
#> Residual standard error: 0.06489 on 183 degrees of freedom
#> Log likelihood =   232, aic =  1268 aicc =  1268, bic(corrected for length) = -5.331
#> 

# As above, another way
spec_a_single  <- x13_spec(spec = "RSA3", transform.function = "Log")
a5 <- x13_automdl(myseries, spec_a_single, identification_estimate.to = "2014-03-01")
a5$regarima
#> y = regression model + arima (1, 0, 2, 0, 1, 1)
#> Log-transformation: yes
#> Coefficients:
#>           Estimate Std. Error
#> Phi(1)     -0.9522      0.030
#> Theta(1)   -1.1092      0.072
#> Theta(2)    0.3625      0.071
#> BTheta(1)  -1.0000      0.045
#> 
#>      Estimate Std. Error
#> Mean  0.02675      0.004
#> 
#> 
#> Residual standard error: 0.06489 on 183 degrees of freedom
#> Log likelihood =   232, aic =  1268 aicc =  1268, bic(corrected for length) = -5.331
#> 

allvar <- pickmdl_data("allvar")

spec_b <- x13_spec(
            spec = "RSA3", transform.function = "Log",
            usrdef.varEnabled = TRUE, 
            usrdef.varType = c("Calendar", "Calendar"), 
            usrdef.var = allvar, 
            outlier.enabled = FALSE, 
            usrdef.outliersEnabled = TRUE,
            usrdef.outliersType = rep("LS", 20), 
            usrdef.outliersDate = c("2009-01-01", "2016-01-01", 
                                    "2020-03-01", "2020-04-01", "2020-05-01", 
                                    "2020-06-01", "2020-07-01", "2020-08-01", 
                                    "2020-09-01", "2020-10-01", "2020-11-01", 
                                    "2020-12-01", "2021-01-01", "2021-02-01",
                                    "2021-03-01", "2021-04-01", "2021-05-01",
                                    "2021-06-01", "2021-07-01", "2021-08-01"))
b <- x13_pickmdl(myseries, spec_b, identification_end = c(2020, 2))                                     
b$regarima
#> y = regression model + arima (0, 1, 1, 0, 1, 1)
#> Log-transformation: yes
#> Coefficients:
#>           Estimate Std. Error
#> Theta(1)   -0.8214      0.045
#> BTheta(1)  -0.8562      0.051
#> 
#>               Estimate Std. Error
#> arb_dag       0.010938      0.001
#> skuddar       0.055092      0.028
#> LS (1-2009)  -0.090294      0.030
#> LS (1-2016)  -0.039840      0.030
#> LS (3-2020)  -0.004550      0.053
#> LS (4-2020)  -0.126323      0.069
#> LS (5-2020)   0.067283      0.068
#> LS (6-2020)  -0.013231      0.068
#> LS (7-2020)   0.108865      0.068
#> LS (8-2020)  -0.118575      0.069
#> LS (9-2020)   0.043492      0.068
#> LS (10-2020) -0.021723      0.068
#> LS (11-2020)  0.027382      0.068
#> LS (12-2020)  0.031743      0.069
#> LS (1-2021)  -0.071634      0.069
#> LS (2-2021)  -0.004503      0.069
#> LS (3-2021)   0.009891      0.069
#> LS (4-2021)   0.024852      0.069
#> LS (5-2021)  -0.051741      0.069
#> LS (6-2021)  -0.035384      0.070
#> LS (7-2021)  -0.017531      0.069
#> LS (8-2021)   0.056976      0.060
#> 
#> 
#> Residual standard error: 0.04928 on 163 degrees of freedom
#> Log likelihood = 290.6, aic =  1180 aicc =  1188, bic(corrected for length) = -5.352
#> 

# automdl instead  
b1 <- x13_automdl(myseries, spec_b, identification_end = c(2020, 2))
b1$regarima
#> y = regression model + arima (0, 1, 1, 0, 1, 1)
#> Log-transformation: yes
#> Coefficients:
#>           Estimate Std. Error
#> Theta(1)   -0.8214      0.045
#> BTheta(1)  -0.8562      0.051
#> 
#>               Estimate Std. Error
#> arb_dag       0.010938      0.001
#> skuddar       0.055092      0.028
#> LS (1-2009)  -0.090294      0.030
#> LS (1-2016)  -0.039840      0.030
#> LS (3-2020)  -0.004550      0.053
#> LS (4-2020)  -0.126323      0.069
#> LS (5-2020)   0.067283      0.068
#> LS (6-2020)  -0.013231      0.068
#> LS (7-2020)   0.108865      0.068
#> LS (8-2020)  -0.118575      0.069
#> LS (9-2020)   0.043492      0.068
#> LS (10-2020) -0.021723      0.068
#> LS (11-2020)  0.027382      0.068
#> LS (12-2020)  0.031743      0.069
#> LS (1-2021)  -0.071634      0.069
#> LS (2-2021)  -0.004503      0.069
#> LS (3-2021)   0.009891      0.069
#> LS (4-2021)   0.024852      0.069
#> LS (5-2021)  -0.051741      0.069
#> LS (6-2021)  -0.035384      0.070
#> LS (7-2021)  -0.017531      0.069
#> LS (8-2021)   0.056976      0.060
#> 
#> 
#> Residual standard error: 0.04928 on 163 degrees of freedom
#> Log likelihood = 290.6, aic =  1180 aicc =  1188, bic(corrected for length) = -5.352
#> 

# effect of identify_t_filter and identify_s_filter
set.seed(1)
rndseries <- ts(rep(1:12, 20) + (1 + (1:240)/20) * runif(240) + 0.5 * c(rep(1, 120), (1:120)^2), 
                frequency = 12, start = c(2000, 1))
spec_c <- x13_spec(outlier.enabled = FALSE)               
c1 <- x13_automdl(rndseries, spec_c, identification_end = c(2009, 12))    
c1$decomposition
#> Monitoring and Quality Assessment Statistics:
#>       M stats
#> M(1)    0.200
#> M(2)    0.024
#> M(3)    0.020
#> M(4)    0.139
#> M(5)    0.268
#> M(6)    0.783
#> M(7)    1.307
#> M(8)    0.915
#> M(9)    0.730
#> M(10)   0.120
#> M(11)   0.120
#> Q       0.472
#> Q-M2    0.535
#> 
#> Final filters: 
#> Seasonal filter:  3x3
#> Trend filter:  13 terms Henderson moving average
c2 <- x13_automdl(rndseries, spec_c, identification_end = c(2009, 12), identify_t_filter = TRUE) 
c2$decomposition
#> Monitoring and Quality Assessment Statistics:
#>       M stats
#> M(1)    0.248
#> M(2)    0.023
#> M(3)    0.000
#> M(4)    0.139
#> M(5)    0.289
#> M(6)    0.768
#> M(7)    1.288
#> M(8)    0.933
#> M(9)    0.731
#> M(10)   0.060
#> M(11)   0.059
#> Q       0.470
#> Q-M2    0.533
#> 
#> Final filters: 
#> Seasonal filter:  3x3
#> Trend filter:  23-Henderson
c3 <- x13_automdl(rndseries, spec_c, identification_end = c(2009, 12), identify_t_filter = TRUE, 
                  identify_s_filter = TRUE)     
c3$decomposition                       
#> Monitoring and Quality Assessment Statistics:
#>       M stats
#> M(1)    0.378
#> M(2)    0.032
#> M(3)    0.069
#> M(4)    0.338
#> M(5)    0.317
#> M(6)    0.751
#> M(7)    1.270
#> M(8)    0.830
#> M(9)    0.735
#> M(10)   0.093
#> M(11)   0.092
#> Q       0.531
#> Q-M2    0.592
#> 
#> Final filters: 
#> Seasonal filter:  3x5
#> Trend filter:  23-Henderson


# Warning when transform.function = "None"
spec_d  <- x13_spec(spec = "RSA3", transform.function = "None")
d <- x13_pickmdl(myseries, spec_d, verbose = TRUE)
#> old_crit2 set to FALSE since no identification specification. 
#> Warning: No model is ok according to criteria
#>  p  d  q bp bd bq 
#>  0  1  1  0  1  1 
#> outlier.from updated: 2021-10-01  No new outliers.

# Warning avoided (when_star) and 2nd (star) model selected 
d2 <- x13_pickmdl(myseries, spec_d, star = 2, when_star = NULL, verbose = TRUE)
#> old_crit2 set to FALSE since no identification specification. 
#>  p  d  q bp bd bq 
#>  0  1  2  0  1  1 
#> outlier.from updated: 2021-10-01  No new outliers.

# automdl since no pickmdl model ok, but still not ok 
d3 <- x13_pickmdl(myseries, spec_d, pickmdl_method = "first_automdl", verbose = TRUE)
#> old_crit2 set to FALSE since no identification specification. 
#> Warning: No model is ok according to criteria
#>  p  d  q bp bd bq 
#>  2  1  1  0  1  1 
#> automdl since no pickmdl model ok
#> outlier.from updated: 2021-10-01  No new outliers.

# airline model (star) since automdl also not ok 
d4 <- x13_pickmdl(myseries, spec_d, pickmdl_method = "first_tryautomdl", verbose = TRUE,
                  when_finalnotok = warning) # also finalnotok warning
#> old_crit2 set to FALSE since no identification specification. 
#> Warning: No model is ok according to criteria
#>  p  d  q bp bd bq 
#>  0  1  1  0  1  1 
#> outlier.from updated: 2021-10-01  No new outliers.
#> Warning: FINAL RUN NOT OK

# As a2, with output = "all"
k <- x13_pickmdl(myseries, spec_b, identification_end = c(2010, 2), output = "all",
                 fastfirst = FALSE) # With TRUE only one model in this case 
k$sa$decomposition  # As a2$decomposition 
#> Monitoring and Quality Assessment Statistics:
#>       M stats
#> M(1)    0.343
#> M(2)    0.349
#> M(3)    1.726
#> M(4)    0.435
#> M(5)    1.870
#> M(6)    0.702
#> M(7)    0.157
#> M(8)    0.197
#> M(9)    0.118
#> M(10)   0.160
#> M(11)   0.153
#> Q       0.610
#> Q-M2    0.646
#> 
#> Final filters: 
#> Seasonal filter:  3x9
#> Trend filter:  23-Henderson
k$mdl_nr            # index of selected model used to identify parameters
#> [1] 1
k$sa_mult[[k$mdl_nr]]$decomposition  # decomposition for model to identify
#> Monitoring and Quality Assessment Statistics:
#>       M stats
#> M(1)    0.264
#> M(2)    0.248
#> M(3)    1.539
#> M(4)    0.000
#> M(5)    3.000
#> M(6)    1.056
#> M(7)    0.285
#> M(8)    0.000
#> M(9)    0.000
#> M(10)   0.000
#> M(11)   0.000
#> Q       0.721
#> Q-M2    0.816
#> 
#> Final filters: 
#> Seasonal filter:  3x5
#> Trend filter:  23-Henderson
k$crit_tab          # Table of criteria 
#>           crit1     crit2      crit3    m_aic
#> [1,] 0.04327429 0.4368229 -0.7117072 344.0338
#> [2,] 0.04325725 0.5369880 -0.7497611 345.7862
#> [3,] 0.04875133 0.5635680  0.0000000 344.3432
#> [4,] 0.04160156 0.9047367 -0.9999670 345.5339
#> [5,] 0.04947505 0.6257011 -0.3132899 346.8342


# Effect of identify_outliers (TRUE is default)
m1 <- x13_pickmdl(myseries, x13_spec("RSA3", outlier.usedefcv = FALSE, outlier.cv = 3), 
                  identification_end = c(2010, 2), identify_outliers = FALSE)
m2 <- x13_pickmdl(myseries, x13_spec("RSA3", outlier.usedefcv = FALSE, outlier.cv = 3), 
                  identification_end = c(2010, 2), identify_outliers = TRUE, 
                  verbose = TRUE, output = "all")
#> old_crit2 set to TRUE based on the date found (2010-03-01) and threshold (2024-10-15).
#>  p  d  q bp bd bq 
#>  0  1  1  0  1  1 
#> outlier.from updated: 2010-03-01  New outliers: 2008-03 
m3 <- x13_pickmdl(myseries, m2$spec, identification_end = c(2018, 2), identify_outliers = TRUE, 
                  verbose = TRUE)
#> old_crit2 set to TRUE based on the date found (2018-03-01) and threshold (2024-10-15).
#>  p  d  q bp bd bq 
#>  2  1  2  0  1  1 
#> outlier.from updated: 2018-03-01  New outliers: 2011-06, 2016-03 

m1$regarima
#> y = regression model + arima (0, 1, 1, 0, 1, 1)
#> Log-transformation: yes
#> Coefficients:
#>           Estimate Std. Error
#> Theta(1)   -0.8434      0.040
#> BTheta(1)  -0.9735      0.044
#> 
#> 
#> Residual standard error: 0.0707 on 185 degrees of freedom
#> Log likelihood =   216, aic =  1285 aicc =  1285, bic(corrected for length) = -5.243
#> 
m2$sa$regarima
#> y = regression model + arima (0, 1, 1, 0, 1, 1)
#> Log-transformation: yes
#> Coefficients:
#>           Estimate Std. Error
#> Theta(1)   -0.8271      0.042
#> BTheta(1)  -0.9189      0.044
#> 
#>             Estimate Std. Error
#> AO (3-2008)  -0.2046      0.068
#> 
#> 
#> Residual standard error: 0.07051 on 184 degrees of freedom
#> Log likelihood = 220.4, aic =  1278 aicc =  1279, bic(corrected for length) = -5.22
#> 
m3$regarima
#> y = regression model + arima (2, 1, 2, 0, 1, 1)
#> Log-transformation: yes
#> Coefficients:
#>           Estimate Std. Error
#> Phi(1)      0.8897      0.145
#> Phi(2)      0.4636      0.080
#> Theta(1)   -0.2235      0.152
#> Theta(2)   -0.2411      0.132
#> BTheta(1)  -0.8802      0.045
#> 
#>             Estimate Std. Error
#> AO (3-2008)  -0.1418      0.054
#> TC (6-2011)  -0.1651      0.040
#> AO (3-2016)  -0.1853      0.054
#> 
#> 
#> Residual standard error: 0.06174 on 179 degrees of freedom
#> Log likelihood = 247.1, aic =  1235 aicc =  1236, bic(corrected for length) = -5.347
#> 


# With corona outliers (even possible when series is not long enough) 
m4 <- x13_pickmdl(myseries, spec_a, verbose = TRUE, corona = TRUE)
#> old_crit2 set to FALSE since no identification specification. 
#>  p  d  q bp bd bq 
#>  2  1  0  0  1  1 
#> outlier.from updated: 2021-10-01  New outliers: 2011-06 
m4$regarima
#> y = regression model + arima (2, 1, 0, 0, 1, 1)
#> Log-transformation: yes
#> Coefficients:
#>           Estimate Std. Error
#> Phi(1)      1.0034      0.062
#> Phi(2)      0.5914      0.062
#> BTheta(1)  -1.0000      0.050
#> 
#>              Estimate Std. Error
#> LS (3-2020)   0.07527      0.065
#> LS (4-2020)  -0.19815      0.092
#> LS (5-2020)   0.05206      0.096
#> LS (6-2020)   0.02156      0.096
#> LS (7-2020)   0.14629      0.100
#> LS (8-2020)  -0.22966      0.102
#> LS (9-2020)   0.15744      0.102
#> LS (10-2020) -0.07074      0.103
#> LS (11-2020)  0.02439      0.103
#> LS (12-2020)  0.09929      0.103
#> LS (1-2021)  -0.19861      0.103
#> LS (2-2021)   0.06577      0.103
#> LS (3-2021)   0.06252      0.103
#> LS (4-2021)  -0.01084      0.103
#> LS (5-2021)  -0.11020      0.103
#> LS (6-2021)   0.04380      0.103
#> LS (7-2021)  -0.04410      0.103
#> LS (8-2021)  -0.03237      0.103
#> LS (9-2021)   0.18665      0.103
#> TC (6-2011)  -0.19161      0.045
#> 
#> 
#> Residual standard error: 0.05883 on 164 degrees of freedom
#> Log likelihood = 248.2, aic =  1263 aicc =  1270, bic(corrected for length) = -5.025
#> 
m5 <- x13_pickmdl(myseries, x13_spec("RSA3", outlier.usedefcv = FALSE, outlier.cv = 3), 
                  identification_end = c(2010, 2), identify_outliers = TRUE, 
                  verbose = TRUE, corona = TRUE) 
#> old_crit2 set to TRUE based on the date found (2010-03-01) and threshold (2024-10-15).
#>  p  d  q bp bd bq 
#>  0  1  1  0  1  1 
#> outlier.from updated: 2010-03-01  New outliers: 2008-03 
m5$regarima 
#> y = regression model + arima (0, 1, 1, 0, 1, 1)
#> Log-transformation: no
#> Coefficients:
#>           Estimate Std. Error
#> Theta(1)   -0.7955      0.048
#> BTheta(1)  -0.9997      0.050
#> 
#>              Estimate Std. Error
#> LS (3-2020)    2.0624      5.852
#> LS (4-2020)  -19.8622      7.464
#> LS (5-2020)    9.9601      7.432
#> LS (6-2020)    2.2136      7.422
#> LS (7-2020)    3.4862      7.422
#> LS (8-2020)  -11.0665      7.432
#> LS (9-2020)   18.5012      7.432
#> LS (10-2020)  -8.2873      7.429
#> LS (11-2020)   2.7988      7.423
#> LS (12-2020)   5.2338      7.423
#> LS (1-2021)  -17.0809      7.423
#> LS (2-2021)    3.4744      7.408
#> LS (3-2021)    4.6064      7.450
#> LS (4-2021)    0.8282      7.464
#> LS (5-2021)  -11.0466      7.432
#> LS (6-2021)    8.4089      7.422
#> LS (7-2021)  -13.3171      7.422
#> LS (8-2021)    4.1312      7.432
#> LS (9-2021)   24.9996      7.432
#> AO (3-2008)  -18.9346      5.543
#> AO (3-2016)  -19.5750      5.543
#> AO (8-2018)   19.7377      5.580
#> AO (4-2017)  -15.8521      5.513
#> TC (10-2018)  14.0022      4.820
#> AO (3-2013)  -19.3536      5.689
#> TC (1-2013)   13.9826      4.881
#> 
#> 
#> Residual standard error: 5.189 on 159 degrees of freedom
#> Log likelihood = -593.7, aic =  1245 aicc =  1256, bic(corrected for length) = 4.073
#> 

 
###########  quarterly series  #############
   
qseries <- pickmdl_data("qseries")    

# Effect of identify_outliers (TRUE is default)
q1 <- x13_pickmdl(qseries, x13_spec("RSA3", outlier.usedefcv = FALSE, outlier.cv = 3), 
                  identification_end = c(2010, 2), identify_outliers = FALSE)
q2 <- x13_pickmdl(qseries, x13_spec("RSA3", outlier.usedefcv = FALSE, outlier.cv = 3), 
                  identification_end = c(2010, 2), identify_outliers = TRUE, 
                  verbose = TRUE, output = "all")
#> old_crit2 set to TRUE based on the date found (2010-07-01) and threshold (2024-10-15).
#>  p  d  q bp bd bq 
#>  0  1  1  0  1  1 
#> outlier.from updated: 2010-07-01  New outliers: 2008-07 
q3 <- x13_pickmdl(qseries, q2$spec, identification_end = c(2018, 2), identify_outliers = TRUE, 
                  verbose = TRUE)
#> old_crit2 set to TRUE based on the date found (2018-07-01) and threshold (2024-10-15).
#>  p  d  q bp bd bq 
#>  0  1  1  0  1  1 
#> outlier.from updated: 2018-07-01  New outliers: 2014-04 

q1$regarima
#> y = regression model + arima (0, 1, 1, 0, 1, 1)
#> Log-transformation: yes
#> Coefficients:
#>           Estimate Std. Error
#> Theta(1)   -0.2941      0.125
#> BTheta(1)  -0.7877      0.104
#> 
#> 
#> Residual standard error: 0.04514 on 59 degrees of freedom
#> Log likelihood = 102.1, aic = 369.5 aicc = 369.9, bic(corrected for length) = -6.063
#> 
q2$sa$regarima
#> y = regression model + arima (0, 1, 1, 0, 1, 1)
#> Log-transformation: yes
#> Coefficients:
#>           Estimate Std. Error
#> Theta(1)   -0.4205      0.120
#> BTheta(1)  -0.7619      0.106
#> 
#>               Estimate Std. Error
#> LS (III-2008)   -0.134      0.037
#> 
#> 
#> Residual standard error: 0.04142 on 58 degrees of freedom
#> Log likelihood = 107.6, aic = 360.6 aicc = 361.3, bic(corrected for length) = -6.168
#> 
q3$regarima
#> y = regression model + arima (0, 1, 1, 0, 1, 1)
#> Log-transformation: yes
#> Coefficients:
#>           Estimate Std. Error
#> Theta(1)   -0.3971      0.122
#> BTheta(1)  -0.7461      0.106
#> 
#>               Estimate Std. Error
#> LS (III-2008) -0.13357      0.035
#> AO (II-2014)   0.09058      0.031
#> 
#> 
#> Residual standard error: 0.03882 on 57 degrees of freedom
#> Log likelihood = 111.7, aic = 354.3 aicc = 355.4, bic(corrected for length) = -6.231
#> 


# With corona outliers (even possible when series is not long enough) 
q4 <- x13_pickmdl(qseries, spec_a, verbose = TRUE, corona = TRUE)
#> old_crit2 set to FALSE since no identification specification. 
#>  p  d  q bp bd bq 
#>  2  1  0  0  1  1 
#> outlier.from updated: 2021-10-01  New outliers: 2008-07 
q4$regarima
#> y = regression model + arima (2, 1, 0, 0, 1, 1)
#> Log-transformation: yes
#> Coefficients:
#>           Estimate Std. Error
#> Phi(1)      0.4350      0.134
#> Phi(2)      0.3402      0.130
#> BTheta(1)  -1.0000      0.103
#> 
#>                 Estimate Std. Error
#> LS (I-2020)   -0.0398499      0.038
#> LS (II-2020)  -0.0287782      0.042
#> LS (III-2020) -0.0745579      0.042
#> LS (IV-2020)   0.0850556      0.043
#> LS (I-2021)   -0.0006202      0.043
#> LS (II-2021)  -0.1151605      0.043
#> LS (III-2021)  0.0442247      0.043
#> LS (III-2008) -0.1253552      0.033
#> 
#> 
#> Residual standard error: 0.0335 on 50 degrees of freedom
#> Log likelihood = 116.8, aic = 358.1 aicc = 364.5, bic(corrected for length) = -6.06
#> 
q5 <- x13_pickmdl(qseries, x13_spec("RSA3", outlier.usedefcv = FALSE, outlier.cv = 3), 
                  identification_end = c(2010, 2), identify_outliers = TRUE, 
                  verbose = TRUE, corona = TRUE) 
#> old_crit2 set to TRUE based on the date found (2010-07-01) and threshold (2024-10-15).
#>  p  d  q bp bd bq 
#>  0  1  1  0  1  1 
#> outlier.from updated: 2010-07-01  New outliers: 2008-07 
q5$regarima 
#> y = regression model + arima (0, 1, 1, 0, 1, 1)
#> Log-transformation: yes
#> Coefficients:
#>           Estimate Std. Error
#> Theta(1)   -0.3110      0.136
#> BTheta(1)  -0.8427      0.084
#> 
#>                Estimate Std. Error
#> LS (I-2020)   -0.030835      0.035
#> LS (II-2020)  -0.021382      0.036
#> LS (III-2020) -0.082770      0.036
#> LS (IV-2020)   0.077653      0.036
#> LS (I-2021)    0.002271      0.037
#> LS (II-2021)  -0.098601      0.037
#> LS (III-2021)  0.034777      0.037
#> LS (III-2008) -0.133177      0.032
#> AO (II-2014)   0.093274      0.027
#> AO (II-2011)   0.085193      0.027
#> 
#> 
#> Residual standard error: 0.03105 on 49 degrees of freedom
#> Log likelihood = 124.8, aic = 344.2 aicc = 351.8, bic(corrected for length) = -6.146
#> 


# Demonstrate strange behavior of x13 with TC at end. Updating outlier.from matters.
# Explains why TC (III-2021) is in q4 and not in q5 
q11 <- x13(window(qseries, end = c(2020, 1)), 
    spec = x13_spec(spec = "RSA3", transform.function = "Log", 
    usrdef.outliersEnabled = TRUE, usrdef.outliersType = "TC", 
    usrdef.outliersDate = "2020-01-01"))  
q12 <- x13(window(qseries, end = c(2020, 1)), # same with outlier.from 
    spec = x13_spec(spec = "RSA3", transform.function = "Log", 
    usrdef.outliersEnabled = TRUE, usrdef.outliersType = "TC", 
    usrdef.outliersDate = "2020-01-01", outlier.from = "2020-04-01")) 
q11$regarima 
#> y = regression model + arima (0, 1, 1, 0, 1, 1)
#> Log-transformation: yes
#> Coefficients:
#>           Estimate Std. Error
#> Theta(1)   -0.2676      0.135
#> BTheta(1)  -0.9100      0.093
#> 
#> 
#> Residual standard error: 0.04206 on 53 degrees of freedom
#> Log likelihood = 94.55, aic = 327.9 aicc = 328.4, bic(corrected for length) = -6.194
#> 
q12$regarima
#> y = regression model + arima (0, 1, 1, 0, 1, 1)
#> Log-transformation: yes
#> Coefficients:
#>           Estimate Std. Error
#> Theta(1)   -0.2693      0.136
#> BTheta(1)  -0.9248      0.092
#> 
#>             Estimate Std. Error
#> TC (I-2020) -0.03069      0.043
#> 
#> 
#> Residual standard error: 0.04166 on 52 degrees of freedom
#> Log likelihood = 94.81, aic = 329.4 aicc = 330.2, bic(corrected for length) = -6.141
#>