Testing the significance of signals

When we work with relatively short and low-frequency time series including linear trends and seasonality, there is the risk of taking the white noise as signal [28]. Therefore, it is useful to test the significance of signals before moving onto further analyses such as DLMs and the ARDL bounds testing. To test the existence of signal against noise, Gershunov et al. [28] proposed a Monte Carlo test approach based on the standard deviations of rolling correlations. In this test, standard deviations of empirical rolling correlations are compared to those of two correlated white noise series replicated many times with the same magnitude of correlation as the empirical series [28]. Higher or lower variability in the standard deviations of rolling correlations than that of two white noise series suggests that the series is not dominated by other characteristics such as trend or seasonality; hence, the correlation between the series is physical. This procedure does not test the existence of a precedence relationship between the series. It rather tests whether the correlation patterns within predefined time windows are significantly different from those in the white noise series. dLagM includes two functions to implement this approach. To the best of our knowledge, there is no readily available software to implement this test other than dLagM.

Rolling correlations are generated and plotted by the following function:

rolCorPlot(x, y, width, start = 1, level = 0.95, main = NULL, SDtest = TRUE)

which implements the Monte Carlo test approach utilizing the roll_cor() function from the roll package [30] and plots the rolling correlations between x and y for the given window lengths by width. It also plots the average rolling correlation over the window lengths after smoothing with nonlinear running median filter by runmed() function of stats package [31] along with level% confidence limits. To test the significance of signal, SDtest is set to TRUE. The following function

sdPercentiles(n = 150, cor = 0.5, width = 5, N = 500, percentiles = c(.05, .95))

produces one-tailed test limits for this approach for the given length n, correlation cor, width and percentiles. The number of Monte Carlo replications is specified by N. If the standard deviations of running correlations between empirical series are outside these limits, we conclude that the signal in the series is significantly more or less variable than expected from noise [28]. sdPercentiles() can be used independently of the rolCorPlot() to produce the limits for any pair of series.

For illustrative purposes, we use the dataset seaLevelTempSOI composed of monthly mean global sea level (GMSL, compared to 1993-2008 average) [32], monthly mean global land-ocean temperature anomalies (1951-1980 as a baseline period) [33], and monthly Southern Oscillation Index (SOI) [34] series between July 1885 and June 2013. We use the last 8 years of the GMSL and temperature anomalies series to demonstrate the usage of the function with the following code chunk.

data(seaLevelTempSOI)

level.ts <- ts(seaLevelTempSOI[1499:1595,1], start = c(2005,5), freq = 12)

temp.ts <- ts(seaLevelTempSOI[1499:1595,2], start = c(2005,5), freq = 12)

rolCorPlot(y = level.ts, x = temp.ts, width = c(3, 5, 7, 11), level = 0.95,

      main = “Rolling correlations between sea levels and temperature”,

      SDtest = TRUE)

The output includes the rolling correlations for each width in $rolCor, the average rolling correlations filtered by running median filter against outliers in $rolcCor.avr.filtered, the unfiltered average rolling correlations in $rolcCor.avr.raw, standard deviations of rolling correlations for each width in rolCor.sd, Pearson correlation between two series in $rawCor, the percentiles of MC distribution of standard deviations of rolling correlations as the test limits in $sdPercentiles, and the standard deviations of rolling correlations for each width along with level% and (1 − level)% limits in $test. The following section of the output shows the test results:

$test

 Width SDrolCor     95%     5%

1   3 0.7431643 0.7486708 0.6505622

2   5 0.5806556 0.5548243 0.4171213

3   7 0.4800409 0.4715809 0.3172123

4  11 0.3975792 0.3733404 0.2047842

The rolling correlation standard deviations (SDs) are outside the limits for the widths 5, 7, and 11. However, the SD for the width of 3 is inside the limits. Thus, the signal between the GMSL and the temperature series is significant for wider window lengths; hence, the rolling correlations between series are not superfluous. Fig 1 shows the plot generated by the rolCorPlot function.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. Rolling correlations between the monthly GMSL and mean land-ocean temperature series.

https://doi.org/10.1371/journal.pone.0228812.g001

In Fig 1, the dashed red lines show the limits of the 95% confidence interval for the mean of the average rolling correlations over the time points, which is shown by the horizontal solid line. The bold, solid line shows the average rolling correlation over the widths. A line is plotted to show the rolling correlations for each width. In accordance with the test results, rolling correlations for widths of 5, 7, and 11 vary more than those for the width of 3. High rolling correlations are recorded between 5 and 7-years-apart sea level and the temperature observations between 2009 and 2010, after both 2010 and 2012 which are significantly different from the rolling correlations between two independent white noise series.

Finite distributed lag models

DLMs are used to model the current and delayed effects of an independent {X_t} series on a dependent {Y_t} series, where i = 1, 2, …, n. The infinite linear DLM is written as follows: where α shows the overall mean effect and β-parameters represent the lag weights. β₀ shows the effect of independent series and β_s shows the effect of lags of the independent series for s = 1, 2, …. The parameter ϵ_t is the stationary error term for the time point t with zero mean and constant variance [35]. When the summation is bounded above by q, the DLM turns out to be a finite linear DLM of order q. Notice that the number of available {Y_t, X_t} pairs for estimation of model parameters is n − q. A weakness of finite DLMs is the existence of multicollinearity with increasing lag orders due to including multiple lags of the same series in the model.

To produce parameter estimates and their significance tests, we employ the dlm() function, which is capable of fitting multiple independent variables with predefined subsets of their lags. Call to the function dlm() is as follows:

dlm(formula,  data,  x,  y,  q,  remove)

The dependent and independent series can be fed into the dlm() function either using formula and data arguments or x and y arguments. If there is only one independent series, use of x and y arguments provides the straightforward way. When there are multiple independent series, all the dependent and independent series should be included in the argument data of the class data.frame and the model should be passed on by the argument formula, which is a usual formula object of R. There is no need to include the lags of independent series in the formula since they are automatically generated by dlm(). The column names of data corresponding to the independent and dependent series must be the same as those in the formula.

The argument q shows the finite lag length. To drop predefined lags of some independent series from the model, we use the remove argument. Each element of the list remove shows particular lags that will be removed from each independent series. To demonstrate the constriction of the remove argument, suppose we arbitrarily have the dependent series Y and independent series X and Z with formula = Y ∼ X + Z and q = 5. Then, if remove is set to NULL, the fitted model will include all lags from 1 to 5 for both X and Z series. To arbitrarily remove the lags 2, 3, and 4 of X and 2 and 4 of Z series, remove should be defined as follows:

> remove <- list(X = c(2,3,4), Z = c(2,4))

Notice that it is possible to fit a model with different lag lengths for each independent series by removing the desired lags using the remove argument and the main series can be removed from the model when 0 is included within the elements of remove.

To illustrate the implementation of the DLMs with dLagM package, we use the GMSL and temperature anomalies series from seaLevelTempSOI dataset. First, we load the dataset to R environment, fit a finite DLM with lag length 7, which gives the minimum BIC, and display the model output with the following code chunk. In accordance with R’s user interface, the model summary, AIC, BIC, residuals, fitted values, and coefficients are displayed by summary(), AIC(), BIC(), residuals(), fitted(), and coef() functions, respectively.

> dlmFit1 <- dlm(formula = GMSL ~ LandOcean,

          data = seaLevelTempOSI, q = 7)

> summary(dlmFit1)

Residuals:

   Min    1Q Median   3Q   Max

-19.8515 -2.5781 0.1053 2.5427 19.9983

Coefficients:

        Estimate Std. Error t value Pr(>|t|)

(Intercept)  -0.008943 0.101086 -0.088 0.9295

LandOcean.t   0.746159 1.153242 0.647  0.5177

LandOcean.1  -0.176851 1.158004 -0.153 0.8786

LandOcean.2  -1.618818 1.154572 -1.402 0.1611

LandOcean.3  -1.977914 1.151842 -1.717 0.0861 .

LandOcean.4  -2.488326 1.151439 -2.161 0.0308 *

LandOcean.5  -2.288856 1.154280 -1.983 0.0475 *

LandOcean.6  -2.248660 1.157732 -1.942 0.0523 .

LandOcean.7  -2.441134 1.153526 -2.116 0.0345 *

---

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 4.028 on 1579 degrees of freedom

Multiple R-squared: 0.01724, Adjusted R-squared: 0.01226

F-statistic: 3.463 on 8 and 1579 DF, p-value: 0.0005765

AIC and BIC values for the model:

    AIC    BIC

1 8942.719 8996.422

Although the model is significant at 5% level of significance (p-value: 0.0005765), it gives a poor fit in terms of the adjusted R-squared of 1.2%. The coefficients for the lags greater than 2 are significant either at 5% or 10% level. The distribution of residuals is nearly symmetric. Then, we remove independent series (with 0 in the remove), its first and second lags, and the intercept (with -1 in the model formula) from the model only to illustrate the use of the remove argument.

> remove = list(LandOcean = c(0,1,2))

> dlmFit2 <- dlm(formula = GMSL ~ -1 + LandOcean,

         data = seaLevelTempSOI, q = 6, remove = remove)

> summary(dlmFit2)

Residuals:

   Min    1Q  Median    3Q   Max

-20.0704 -2.6089  0.0976 2.5265 20.1688

Coefficients:

       Estimate Std. Error t value Pr(>|t|)

LandOcean.3  -2.072  1.109 -1.868  0.0619 .

LandOcean.4  -2.266  1.116 -2.030  0.0425 *

LandOcean.5  -2.147  1.116 -1.923  0.0547 .

LandOcean.6  -2.480  1.108 -2.237  0.0254 *

---

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 4.034 on 1585 degrees of freedom

Multiple R-squared: 0.01345, Adjusted R-squared: 0.01096

F-statistic: 5.4 on 4 and 1585 DF, p-value: 0.0002555

AIC and BIC values for the model:

    AIC    BIC

1 8947.933 8974.788

With this model, we get similar results on the significance of the coefficients and the reduced model. Note that since R uses a different formula for R-squared when there is no intercept in the model, R-squared values are not comparable. Also, these outputs are given neither to compare two models nor to draw decisive inferences on the relationship between the temperature anomalies and the sea level series. The intention here is to illustrate the usage of the functions. The residuals, fitted values, and coefficients of the model dlmFit2 are obtained by residuals(dlmFit2), fitted(dlmFit2), and coef(dlmFit2), respectively.

It is possible to impose a polynomial shape on the lag weights to fit a polynomial DLM that is fitted by the following call to the function polyDlm():

polyDlm(x, y, q, k, show.beta = TRUE)

Currently, there is only one independent variable allowed in polyDlm() function. The arguments x and y are used to send independent and dependent series into the function, respectively. The arguments q and k are used to specify the lag and polynomial orders, respectively. To show the values of actual β parameters and the results of their significance tests, show.beta needs to be set to TRUE.

Koyck distributed lag models

Another approach to specify lag weight is to impose a geometrically decreasing pattern on the lag weights. The following equation shows the general form of a geometric DLM with geometrically decreasing lag weights β_s = βϕ^s:

The implementation of geometric DLMs is not practical due to the infinite number of parameters and nonlinearity of the model. Koyck transformation is an efficient way of dealing with these drawbacks of the geometric DLMs. The Koyck form of the geometric DLM is as follows [35]:

As the result, we get a finite DLM after the Koyck transformation. Notice that the first lag of the dependent series Y_t−1 and the error term ν_t are correlated with each other and the error term ν_t depends on both ϵ_t and ϵ_t−1. Therefore, we use instrumental variables approach to estimate model parameters in dLagM [36]. The function koyckDlm() implements the instrumental variables estimation. The call to the function is as follows:

koyckDlm(x, y)

Currently, only one independent variable is allowed in koyckDlm() function. The arguments x and y are used to send the independent and dependent series into the function, respectively. Diagnostics for the Koyck DLM are displayed by

summary(KoyckFit, diagnostics=TRUE)

given the fitted model is KoyckFit. All the diagnostic tests are implemented by the package AER [36]. In the Weak Instruments line of the output, an F test of the first stage regression for weak instruments is displayed. In the second line, the Wu-Hausman test is displayed for the significance of the correlation between an explanatory variable and the error term (endogeneity). In the Koyck DLM, the significance of the correlation between the error term and the lagged dependent variable is not clear. The Wu-Hausman test, which compares the consistency of instrumental variable estimates to ordinary least squares estimates, provides a way for testing the significance of the correlation between the error term and the dependent series. If the endogeneity is present, the instrumental variable estimator is concluded to be consistent while the ordinary least squares estimator is inconsistent. But this test should be used cautiously with small samples due to its asymptotic validity [37].

We fit the Koyck model with GMSL and land-ocean temperature anomalies series and display the summary statistics with the following code chunk:

> KoyckFit <- koyckDlm(x = seaLevelTempSOI$LandOcean, y = seaLevelTempSOI$GMSL)

> summary(KoyckFit, diagnostics=TRUE)

Residuals:

    Min     1Q   Median     3Q     Max

-16.74923 -1.75048 -0.02326 1.68498  11.97262

Coefficients:

       Estimate Std. Error t value Pr(>|t|)

(Intercept) 0.008036 0.070254 0.114 0.909

Y.1      0.733753 0.017582 41.734 <2e-16 ***

X.t     -2.553163 3.607947 -0.708 0.479

---

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.805 on 1591 degrees of freedom

Multiple R-Squared: 0.5297, Adjusted R-squared: 0.5291

Wald test: 906.6 on 2 and 1591 DF, p-value: < 2.2e-16

Diagnostic tests:

         df1  df2    statistic    p-value

Weak instruments 1 1591  70.3332073 1.088953e-16

Wu-Hausman     1 1590   0.9925055 3.192824e-01

                 alpha    beta    phi

Geometric coefficients: 0.03018417 -2.553163 0.7337525

We have a significant model and a coefficient with lower adjusted R-squared. There is no significant endogeneity by the Wu-Hausman test and no significant weak instruments problem by the Weak Instruments test.

Autoregressive distributed lag models

The autoregressive DLMs are employed to implement a regression analysis between a dependent series and k number of independent series. When there is only one independent series, an ARDL model of orders p and q is denoted by ARDL(p, q), which consists of p lags of independent and q lags of dependent series. The lags of the dependent series make the model autoregressive. When the number of lags of ith independent series is shown by p_i, i = 1, …, k, we write the ARDL model as in Eq (1): (1) where μ₀ is the constant, Y_t and X_ti are respectively dependent and ith independent series, p_i is the lag order of ith independent series, q is the autoregressive order of the model, and e_t shows the innovations.

To implement autoregressive DLMs, the call to the ardlDlm() function is

ardlDlm(formula = NULL, data = NULL, x = NULL, y = NULL, p = 1, q = 1, remove = NULL)

The same convention as the dlm() function is followed to feed the data and the model into the ardlDlm() function. The argument p shows the values of lag orders such that p₁ = p₂ = ⋯ = p_k = p. The argument q shows the value of autoregressive orders of the model. This makes the call to the function easier when all lag orders are equal. When the lag order of each independent series differs, we employ remove argument to remove unwanted lags. It is also used to remove some lags of dependent series, simultaneously. In this sense, it is very flexible in fitting parsimonious models. It has the following list structure:

remove <- list(p = list(), q = c())

Each element of the list p shows particular lags that will be removed from each independent series. To remove the main series from the model or to fit an ARDL(0,q) model, we include 0 within the elements of p. The element q is a vector showing the autoregressive lags of dependent series to be removed. Following the example given for DLM implementation, to remove the main series of X and Z, the second lag of X and the first autoregressive lag of Y from the model, we define remove as follows:

remove <- list(p = list(X = c(0,2), Z = c(0)), q = c(1)).

Note that the generic functions mentioned for dlm() function works for all the models fitted by dLagM.

To demonstrate the use of the function with multiple independent series, we use all three series available by the seaLevelTempSOI dataset available in the package. The following code chunk fits an ARDL model for the GMSL (Y) series with land-ocean temperature anomalies (X₁) and SOI (X₂) series. For demonstration purposes, we arbitrarily take p₁ = 2, p₂ = 1, q = 4 and fit the ARDL model.

> formula1 <- GMSL ~ LandOcean + SOI

> remove <- list(p = list(SOI = c(2)))

> ARDLfit1 <- ardlDlm(formula = formula1, data = seaLevelTempSOI, p = 2, q = 4, remove = remove)

> summary(ARDLfit1)

Residuals:

   Min    1Q Median   3Q   Max

-11.0820 -1.2365 -0.0346 1.1360 7.5509

Coefficients:

        Estimate Std. Error t value Pr(>|t|)

(Intercept)   0.003920 0.049067 0.080 0.93633

LandOcean.t   1.129695 0.533647 2.117 0.03442 *

LandOcean.1   0.291991 0.545341 0.535 0.59243

LandOcean.2  -1.369952 0.531626 -2.577 0.01006 *

SOI.t      -0.010608 0.006061 -1.750 0.08027 .

SOI.1      0.002837 0.006078  0.467 0.64071

GMSL.1      0.951282 0.024389 39.004 < 2e-16 ***

GMSL.2      0.081073 0.028575  2.837 0.00461 **

GMSL.3     -0.745129 0.028568 -26.083 < 2e-16 ***

GMSL.4      0.240132 0.024381  9.849 < 2e-16 ***

---

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.957 on 1581 degrees of freedom

Multiple R-squared: 0.7718,Adjusted R-squared: 0.7705

F-statistic: 594 on 9 and 1581 DF, p-value: < 2.2e-16

The coefficient related to X₁ series and its second lag, and all four lags of the Y series are significant at the 5% significance level. Also, X₂ series is significant at the 10% level. The model is significant the 5% level with a p-value less than 2.2e-16. We can take a further step and remove insignificant coefficients from the model as follows:

> remove <- list(p = list(LandOcean = c(1), SOI = c(1,2)))

> ARDLfit2 <- ardlDlm(formula = formula1, data = seaLevelTempSOI, p = 2, q = 4, remove = remove)

> summary(ARDLfit2)

Residuals:

    Min   1Q Median  3Q  Max

 -11.060 -1.237 -0.030 1.159 7.504

Coefficients:

       Estimate Std. Error t value Pr(>|t|)

(Intercept)   0.003949 0.049043 0.081 0.93583

LandOcean.t   1.185116 0.521855 2.271 0.02328 *

LandOcean.2  -1.320628 0.519258 -2.543 0.01108 *

SOI.t      -0.008862 0.004722 -1.877 0.06074 .

GMSL.1      0.951720 0.024314 39.142 < 2e-16 ***

GMSL.2      0.080261 0.028509  2.815 0.00493 **

GMSL.3     -0.744486 0.028536 -26.089 < 2e-16 ***

GMSL.4      0.239964 0.024349  9.855 < 2e-16 ***

---

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.956 on 1583 degrees of freedom

Multiple R-squared: 0.7717, Adjusted R-squared: 0.7707

F-statistic: 764.4 on 7 and 1583 DF, p-value: < 2.2e-16

According to this model, the previous year’s Southern Oscillation Index, temperature anomaly recorded in the current year and two years before, and sea levels of one, two, three, and four years before have a significant impact on the sea level in a given year. The model is significant at the 5% level (P < 2.2e − 16), the adjusted R-squared value is at a moderate to high level with 77.07%. A model selection procedure needs to be applied to find the optimum model for this data. Here, we aim to explain the use of the arguments of the ardlDlm() function.

Forecasting with distributed lag models

For all the models implemented by the dLagM package, forecasts are generated by the generic forecast() function:

forecast(model, x, h = 1, interval = FALSE, level = 0.95, nSim = 500)

The argument model is a fitted model by one of the functions of dLagM package. x includes the new data for independent series with a forecast horizon of h. The (1 − level)% prediction intervals for forecasts are obtained by the Monte Carlo approach using a Gaussian error distribution with zero mean and empirical variance of the dependent series if interval is set to TRUE.

To illustrate the use of the function, the following code chunk generates 3 months ahead sea level forecasts with 95% prediction intervals using the fitted ARDL model ARDLfit2.

> x <- matrix(c(0.06, 0.07, 0.08, 13.9, 8.1, -0.5), nrow = 2, ncol = 3, byrow = T)

> forecast(model = ARDLfit2, x = x, h = 3, interval = TRUE)$forecasts

   95% LB  Forecast  95% UB

1 -1.565833  2.371956 5.948790

2 2.414651  6.052777 9.696141

3 -7.784304 -2.816127 2.278009

Here, the matrix x provides the new values of adjusted temperature anomalies and SOI values in the forecast horizon of 3. When there are multiple independent series, x must be entered as a matrix, each row of which has h elements and corresponds to each independent series.

Goodness-of-fit measures

There are three auxiliary functions included in dLagM: finiteDLMauto(), GoF(), and sortScore(). The GoF() function computes mean absolute error (MAE), mean squared error (MSE), mean percentage error (MPE), symmetric mean absolute percentage error (sMAPE), mean absolute percentage error (MAPE), mean absolute scaled error (MASE), mean relative absolute error (MRAE), geometric mean relative absolute error (GMRAE), mean bounded relative absolute error (MBRAE), unscaled MBRAE (UMBRAE) for distributed lag models using the formulations of Hyndman et al. [38, p. 26] and Chen et al. [29]. When there is only one model passed in the function, MBRAE and UMBRAE are not calculated. The function sortScore displays sorted AIC, BIC, MASE, MAPE, sMAPE, MRAE, GMRAE, or MBRAE scores in the ascending order to flag the best performing model easily.

The following code chunk computes the sorted MASE and other measures over the fitted DLMs in the previous sections.

> sortScore(x = MASE(dlmFit1, dlmFit2, KoyckFit, ARDLfit2), score = “mase”)

        n    MASE

ARDLfit2 1591 0.6561139

KoyckFit 1594 0.9399988

dlmFit2  1589 1.3694767

dlmFit1  1588 1.3702182

> GoF(dlmFit1, dlmFit2, KoyckFit, ARDLfit2)

       n    MAE    MPE   MAPE   sMAPE    MASE

dlmFit1 1588 3.131399 0.8207351 1.353352 1.6745267 1.3702182

dlmFit2 1589 3.133837 0.8035919 1.315302 1.6875275 1.3694767

KoyckFit 1594 2.161233 0.5037362 3.065148 1.0002000 0.9399988

ARDLfit2 1591 1.502217 0.3403850 2.367279 0.7693254 0.6561139

         MSE   MRAE   GMRAE   MBRAE   UMBRAE

dlmFit1 16.13457 7.890835 1.4464045 1.8943210 -2.1181668

dlmFit2 16.23229 8.740936 1.4243370 0.5903609 1.4411730

KoyckFit 7.85227 3.864511 0.9760243 0.4545904 0.8334846

ARDLfit2 3.80717 7.013070 0.6688625 0.1175530 0.1332126

Specification of finite lag length

Specification of lag length for finite DLMs is an issue that needs to be elaborated. dLagM includes finiteDLMauto() function to help finding the best model according to goodness-of-fit (GOF) statistics. It automatically searches for the optimal model in terms of one of MASE, AIC, BIC, GMRAE, MBRAE, or adjusted R-square. It is capable of handling multiple independent series for searching over finite DLMs and it works with only one independent series for polynomial DLMs. The call to the function is as follows:

finiteDLMauto(formula, data, x, y, q.min = 1, q.max = 10,

        k.order = NULL, model.type = c(“dlm”,“poly”),

        error.type = c(“MASE”,“AIC”,“BIC”,“GMRAE”,

        “MBRAE”, “radj”), trace = FALSE)

The input structures for formula, data, x, and y arguments are the same as those for the dlm() function. The arguments q.min and q.max specify the bounds of lag order q for the search. When the model.type is set to poly, the arguments x and y must be used to send the independent and dependent series into the function. The function fits all finite linear or polynomial DLMs for all lag and/or polynomial orders and calculates or extract MASE, AIC, BIC, GMRAE, MBRAE, adjusted r-squared, and the p-value of Ljung-Box test from the model objects. Then, it sorts the models according to the desired measure. By this way, it displays the effect of lag and/or polynomial order on the error/GOF measures to help the user to determine the optimal value of orders of the model.

Specification of orders for ARDL bounds test

Selection of orders for ARDL bounds test has a direct impact on the resulting inferences. The number of possible orders quickly increases as the number of independent series (k), p_i and q increase. For instance, for k = 3 and p₁ = p₂ = p₃ = q = 5, the total number of possible models is 1080. In some economic applications, the orders are determined by the theory. However, in other cases, the mainstream approach to find a suitable set of orders is based on searching over the full-set of possible models using AIC, BIC, HQIC (Hannan-Quinn IC), or adjusted r-square statistics. We call this approach “full-search” in the rest of this section. This approach is prone to the weaknesses of the used statistic, and its computer implementation is not quite time-efficient and troublesome for practitioners with limited programming skills. The mainstream econometric software EViews includes an add-in program to run the full-search approach [27, 47]. In dLagM package, in addition to the full-search, we implement a two-staged approach to the selection of the orders by the ardlBoundOrders() function. This function provides the users at all levels of software skills with a user-friendly tool for the specification of orders. The main advantages of the two-staged approach are that it utilizes MASE, and GMRAE in addition to AIC and BIC as the decision criterion for the search and it is more time-efficient when the values of k, p_i, and q are large. This makes it possible to search with a larger maximum limit on the orders for long series. For advanced users, the proposed algorithm can be coded under any R package or software that fits ARDL models as an add-in program or package such as EViews.

The ardlBoundOrders() function computes the optimal lag structure for the short-run relationships and autoregressive part of the ARDL model for the ARDL bounds testing. The call to the function is as follows:

ardlBoundOrders(data = NULL, formula = NULL, ic = c(“AIC”, “BIC”,

         “MASE”, “GMRAE”), max.p = 15, max.q = 15,

         FullSearch = FALSE)

When FullSearch = TRUE, the full-search is implemented, and when it is set to FALSE, the two-staged search approach is implemented as in Algorithm 1.

Algorithm 1.

1. Assume all p-orders are equal (p₁ = p₂ = ⋯ = p_k = p), set p = 1 and q = 1.
2. Fit the conditional error correction model. Increase the value of q by one. Repeat this step until max.q is reached. Then, move on to the next step.
3. If p is less than or equal to max.p, increase the value of p by one and go to step 2. Otherwise, go to the next step.
4. Find the (p*, q*) pair that gives the optimal value to the stat.
5. Expand a grid with all combinations of k factors with p* levels.
6. Take each combination in the grid as the value of p₁, p₂, …, p_k and fit the conditional error correction model with orders (p₁, p₂, …, p_k, q*).
7. Report the combination that optimizes the stat as the optimal value or orders for the ARDL bounds test.

The first stage of Algorithm 1 is composed of the steps 1-3, where we find the number of autoregressive lags (q) that gives us a suitable value for the desired GOF statistic when all the orders of short-run relationships are assumed to be equal. After fixing the value of q at the first stage, we find the optimal combination of the p-orders at the second stage. To benchmark the computational efficiency of our approach with the full-search, suppose we have k = 3, max.p = 10 and max.q = 10 which define a moderate-sized search space. With the full-search approach, we need to fit 11³ ⋅ 10 = 13, 310 models. Whereas in our approach, we fit 11 ⋅ 10 = 110 models at the first stage. The number of models to be fitted at the second stage depends on the optimal value of p. In the worst-case scenario, we get p = max.p and fit 11³ = 1, 331 models at the second stage. In total, the maximum number of models to be fitted in our approach is 1,441 which is nearly 10% of the full-search approach. Therefore, in any platform it is implemented, our approach will run faster than the full-search.

The main drawback of our approach is the loss of information because of not searching the whole space of possible models. However, since we find an optimal value for q by assuming all p-orders are equal, this takes us to the proximity of the optimal set of lag orders in our search algorithm. Then, at the second stage, we fine-tune by allowing the p-orders to vary across the independent series in the short-run relationship part of the ARDL model. The full-search algorithm gives the optimal values of each GOF statistic as well as the corresponding lag orders to benchmark with the same values resulting from our algorithm to see the magnitude of information loss. To demonstrate this in practice, we run both full-search and proposed search approaches on 4 yearly, quarterly, monthly, and daily datasets and compare the values of GOF statistics and lag orders with k = 2, 3, and various settings of max.p and max.q. The following benchmarking code is implemented on a MacBook Pro computer with 2.9 GHz Intel Core i7 processor and 16 GB 2133 MHz LPDDR3 memory.

tic(“Standard”)

orders1<- ardlBoundOrders(data = seaLevelTempSOI, formula = formula1, ic = “AIC”, max.p = 10, max.q = 10, FullSearch = TRUE)

toc()

orders1$p; orders1$q; orders1$min.Stat

tic(“Proposed”)

orders2 <- ardlBoundOrders(data = seaLevelTempSOI, formula = formula1, ic = “AIC”, max.p = 10, max.q = 10)

toc()

orders2$p; orders2$q; orders2$min.Stat

We run this code chunk with AIC, BIC, MASE, and GMRAE statistics and record implementation times, optimal orders, and the value of the GOF statistics for the datasets and model formulas presented in Table 1. The results of trials with two independent series are reported in Table 2 and those with three independent series are presented in Table 3.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Description of the test cases for the comparison of the full-search and the proposed algorithm for the specification of orders in the ARDL bounds testing.

https://doi.org/10.1371/journal.pone.0228812.t001

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 2. Benchmarking of p and q lag orders, implementation time, and the value of GOF statistic for the full-search and the proposed algorithm for Trials 1, 2, 4, and 5.

https://doi.org/10.1371/journal.pone.0228812.t002

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 3. Benchmarking of p and q lag orders, implementation time, and the value of GOF statistic for the full-search and the proposed algorithm for Trial 3.

https://doi.org/10.1371/journal.pone.0228812.t003

For all GOF measures, as expected, our algorithm is 3 to 11 times faster than the full-search algorithm. When AIC is used as the target GOF measure, we get the same or more parsimonious orders with near-optimal AIC values via our search algorithm for all trials. For BIC, we find the same orders in all trials faster than the full-search algorithm. Targeting for MASE produces different lag orders in Trials 1 and 3. In Trial 1, the proposed algorithm recommends a similar but more parsimonious model with a very close MASE to its full-search counterpart. Similarly, the proposed algorithm identifies the same orders as the full-search approach quicker in Trials 1 and 5 with GMRAE. In Trials, 2, 4, and 5, minimizing GMRAE returns more parsimonious models with the proposed algorithm. Overall, this benchmarking shows us that the information loss with the two-stage algorithm is negligible while the gain in computational time is highly significant.

The full-search approach gives optimal orders in terms of the GOF measures. However, optimising the GOF measures and simultaneously keeping the model as parsimonious as possible is challenging and can be done by a trial-and-error method if the user is experienced. The values of GOF measures obtained by our search algorithm is very close to the optimal values by the full-search and the proposed orders constitute more parsimonious models in a computationally efficient way. Another consideration is that the trial-and-error method cannot readily be embedded in artificial intelligence algorithms such as recommender systems. In this sense, our algorithm is useful and efficient when it is embedded in self-learning systems to specify the orders automatically. Also, it provides less experienced users with an efficient tool for order specification.

The dLagM package can be compared to EViews and dynamac packages in terms of its functionality for ARDL bounds testing. We summarize this comparison as follows:

• The main distinction of dLagM is that it implements both the full-search and the proposed algorithms to optimize AIC, BIC, MASE, or GMRAE, and reports the optimized orders for ARDL bounds testing. Both EViews and dynamac do not offer any built-in functions for this aim by default. There is an add-in program available for Eviews to implement the full-search which is limited to a maximum of 10 lags [47].
• dLagM readily displays Breusch-Godfrey, Ljung-Box, Breusch-Pagan, Shapiro-Wilk’s, and Ramsey’s RESET tests as residual diagnostics, and recursive CUSUM, recursive CUSUM of Squares, and recursive MOSUM plots as stability diagnostics. EViews has Jarque-Bera, Breusch-Pagan-Godfrey Serial Correlation LM, Harvey, Glejser, ARCH, and White heteroskedasticity tests as residual diagnostics, and recursive CUSUM, recursive CUSUM of squares, recursive coefficients, N-step forecast tests, and Ramsey’s RESET tests as stability diagnostics. dynamac implements Breusch-Godfrey, Chi-square, and F tests can be implemented by a separate function call. The dynamac::pssbounds() function only displays the F-statistic and critical values for the test.
• As input dLagM needs the lag orders and dataset to run the ARDL bounds test. This is similar in EViews. Whereas dynamac requires fitting the right ARDL model outside the dynamac::pssbounds() function.

Overall, as open-source software, dLagM has almost a similar functionality for ARDL bounds testing as EViews plus some different diagnostic tests and an alternative approach to the specification of orders. It is very easy to hook up the dLagM to other research software in R or other software such as Matlab and Python that can communicate with R due to dLagM’s input/output interface and R’s functionality.

Illustration of ARDL bounds testing with dLagM package

To illustrate the use of the package for ARDL bounds testing, we will investigate short and long-term relationships between change in the sea level (relative to 1993-2008 average), land and ocean temperature anomalies and SOI provided by the seaLevelTempSOI dataset.

In practice, we have to check if the signals in both series are significant or white noise by using the rolCorPlot() function and we need to check the existence of a spurious correlation between the dependent and independent series. Since our aim here is to illustrate the use of the ardlBound() function, we do not present these preliminary analyses here. They are included in the R scripts given in the “Supporting information” section. We specify lag orders using ardlBoundOrders() function as follows:

> orders2 <- ardlBoundOrders(data = seaLevelTempSOI, formula = formula1, ic = “BIC”, max.p = 10, max.q = 10)

> orders2$p

 LandOcean SOI

1      1 1

> orders2$q

[1] 10

From this analysis, we get p₁ = p₂ = 1, q = 10. We put these values into a data.frame object to pass on the main function. Some software starts the summations in Eq (2) from 0. However, since we start them from 1 following Pesaran et al. [12], we add 1 to the elements of the data.frame object:

> p <- data.frame(orders2$q, orders2$p) + 1

Then, we run the ARDL bounds test by the following call:

> test1 <- ardlBound(data = seaLevelTempSOI, formula = formula1, case = 1, p = p, ECM = TRUE)

Here, we fit the model with no intercept and no trend represented by case = 1. Since the full output of this call is too long, we leave running the codes given in “Supporting information” section and getting the full output to readers. The output shows us the Breusch-Godfrey, Ljung-Box, Breusch-Pagan, Shapiro-Wilk’s, and Ramsey’s RESET tests as residual diagnostics along with the results of the ARDL bound test. We get non-normal residuals and significant autocorrelations in the residuals by the Shapiro-Wilk and Breusch-Godfrey tests, respectively. RESET test indicates that there is no model specification issue. Since we find an F-statistic of 118.67, which is greater than all the critical values for all α-levels, we conclude the significance of cointegration between the change sea level, land and ocean temperature anomalies, and SOI. Then, the output shows the coefficient estimates and significance tests for the error correction model. The error correction coefficient is negative (-0.9412) as expected and highly significant (p-value < 2e-16), which is consistent with the ARDL bounds test. This means that there is a significant long-run relationship or cointegration between the change in sea levels and the exploratory series. The last bit of the output shows the values of the long-run coefficients for the land-ocean temperature anomalies and SOI. To get the significance tests for the long-run coefficients, we use summary(test1$model$modelFull).

The stability of the model is another important issue to check in ARDL bounds testing. If both ECM = TRUE and stability = TRUE, the cumulative sum (CUSUM), CUSUM of squares, and MOSUM charts [48, p. 213] [40] are generated using recursive residuals. While the recursive CUSUM, recursive CUSUM of squares plots are also generated by EViews, dLagM package provides recursive MOSUM chart to give more insight into the investigation of model stability issue in ARDL bounds testing. Recursive CUSUM is produced over all previous observations, whereas recursive MOSUM uses a rolling window of previous observations [48]. So, MOSUM is more sensitive to local departures from stability. Recursive MOSUM chart requires a relatively long series to be displayed. Recursive CUSUM of squares can be used to identify instability in the variances of coefficients over time. If there is a structural change in the volatility of the model coefficients, we expect the recursive CUSUM of squares cross the limits. Fig 2 displays the recursive CUSUM, recursive CUSUM of squares, and recursive MOSUM plots for the example.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Recursive CUSUM, CUSUM of squares and MOSUM plots for the model stability.

https://doi.org/10.1371/journal.pone.0228812.g002

Since the recursive CUSUM and MOSUM residuals remain close to 0 without fluctuating outside the 5% limits, there is no sign of change in the coefficients over time. However, recursive CUSUM of squares test indicates that there is an instability in the volatility of coefficients as recursive CUSUM of squares goes outside the 5% limits. This is a good example of instability in the volatility of coefficients while they are stable over time.

Overall, for this example, the model specification should be changed to get acceptable diagnostic results before proceeding with the decision provided by the ARDL bounds test. Considering that the orders used for the test are minimizing BIC among all possible models with max.p = 10 and max.q = 10, minimizing the GOF measures to specify the orders does not necessarily ensure getting stable results. Therefore, it is crucial to check residual and stability diagnostics.