A Forecasting Evaluation Library for
J Demetra+
David de Antonio Liedo∗
January 26, 2017
Abstract
This note presents the main tools currently available in JDemetra+ for the
evaluation of forecasting errors. We explain the main statistical concepts and their
implementation.
Key words: JDemetra+Nowcasting, forecasting evaluation, dynamic factor mod-
els, Diebold-Mariano Test, fixed-smoothing asymptotics.
Contents
Evaluating Forecasting Accuracy: Concepts 2
Diebold-Mariano Test . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
Encompassing Test . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
Efficiency: Bias Test . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
Efficiency: Autocorrelation Test . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
Implementation in J Demetra+ 6
Structure of the library . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
A simple example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
Example in the context of nowcasting . . . . . . . . . . . . . . . . . . . . . . . . 9
∗
National Bank of Belgium, R&D Statistics, e-mail: [Link]@[Link].
1
Evaluating Forecasting Accuracy: Concepts
The prediction errors are defined with a reference i to the information set available at
the time the forecast was made: et|i = yt − ŷt|Fi , where Fi may include lags of yt and
many other variables that do not necessarily refer to the same period. In practice, the
information that will be actually used may be a small subset of Fi .
The properties of these forecast errors can be assessed in isolation or relative to a
benchmark, which we will define as ĕt|i = yt − y̆t|Fi . The benchmark may be a naive
forecast, e.g. random walk, in which case y̆t|Fi would be equal to y̆t|yt−1 = yt−1 . However,
the benchmark could also be a prediction that is regularly published by a forecasting
institute or market analysts, which is not necessarily model-based. In that case, y̆t|Fi
would be given by methods and a subset of Fi which is unknown to us.
For model-based forecasts, we use the following notation: ŷt|Fi = Eθ [yt |Fi ] to highlight
the fact that they are based on model-consistent expectations given by the parameter
vector θ.
In forecasting comparisons involving competing forecasts that result from the same
information set, the subindex i will be removed because it does not play a role. We will
first test the following null hypotheses involving forecast errors:
U nbiasedness : E[et ] = 0 (1)
Autocorrelation : E[et et−1 ] = 0 (2)
Equality in squared errors : E[e2t − ĕ2t ] = 0 (3)
Equality in absolute errors : E[|et | − |ĕt |] = 0 (4)
F orecast ŷt encompasses y̆t : E[(et − ĕt )et ] = 0 (5)
F orecast y̆t encompasses ŷt : E[(ĕt − et )ĕt ] = 0 (6)
The null hypothesis may be rejected in favour of either one or two-sided alternatives.
For example, the hypothesis of equality in forecast errors originally proposed by Diebold
and Mariano (see below) was tested against the two-sided alternative1 (i.e. 6= 0).
An overview of the tests can also be found in Table A.1.
1
In turn, one may decide to define a rejection only when the loss differential e2t − ĕ2t is negative, i.e.
our squared forecast errors are smaller than those coming from the benchmark.
2
Diebold-Mariano Test
The test originally proposed by Diebold and Mariano (1995) considers a sample path of
loss differentials {dt }Tt=1 . In the case of a squared loss function, we have dt = e2t − ĕ2t .
Under the assumption that the loss differential is a covariance stationary series, the sample
¯ converges asymptotically to a normal distribution:
average, d,
√
T d¯ →
d N (µ, 2πfd (0))
− (7)
In particular, they proposed to test the null hypothesis that the forecast errors coming
from the two forecasts bring about the same loss: E[e2t − ĕ2t ] = 0 against the two-sided
alternative. Thus, the resulting p-values represent the probability of obtaining the realized
forecast error differential or a more extreme one in a new experiment if the null hypothesis
was actually true. The test-statistic that will be used to calculate our p-values is computed
as follows:
d¯
DM = s (8)
2π fˆd (0)
T
¯ Consider 2π fˆd (0) = P(T −1)
where 2π fˆd (0) is a consistent estimate of the variance of d. τ =−(T −1) wτ γd (τ ),
1 T
P ¯ t−|τ | − d).
¯ Under the assumption that γd (τ ) = 0 for
where γd (τ ) = (dt − d)(d
T t=|τ |+1
τ ≥ h, we can use a rectangular lag window estimator by setting wτ = 0 for τ ≥ h.
Another option is to use the Heteroscedasticity and Autocorrelation Consistent (HAC)
estimator proposed by Newey and West (1987). In this case, the weights could be given
τ
by a triangular window, wτ = 1 − for τ < h. In this case, however, the consistency
h
property only remains valid when the truncation lag h or bandwidth is a function of the
sample size T .
The idea is to test the statistical significance of the regression of e2t −ĕ2t on an intercept.
In order to determine the statistical significance of the intercept, its associated standard
errors need to take into account the autocorrelation patterns of the regression error, which
are considered in the denominator of equation (8). JDemetra+ exploits the same unified
framework to conduct all tests listed in Table A.1. But given the small sample sizes that
are typical in real-time forecasting applications, which leads to an over-rejection of the
3
Table A.1: Forecasting Evaluation Tests
Test Null hypothesis Statistic Asym. Finite JDemetra+
theory sample implementation
ē
Bias E[et ] = 0 B= s N (0, 1) KV(2005) BiasTest
2π fˆe (0)
T
ρ̄
Autocorrelation E[et et−1 ] = 0 AR = s N (0, 1) KV(2005) EfficiencyTest
2π fˆρ (0)
T
4
d¯
Diebold-Mariano dt ≡ L1,t − L2,t = 0 DM = s N (0, 1) KV(2005) DieboldMarianoTest
2π fˆd (0)
T
d¯1
Encompassing 1 de1,t ≡ E[(et − ĕt )et ] = 0 E1 = s N (0, 1) KV(2005) EncompassingTest
2π fˆd1 (0)
T
d¯2
Encompassing 2 de2,t ≡ E[(ĕt − et )ĕt ] = 0 E2 = s N (0, 1) KV(2005) EncompassingTest
2π fˆd2 (0)
T
null hypothesis, we follow the fixed-smoothing (FS) asympotics proposed by Coroneo and
Iacone (2015) exploiting the finite sample distributions of Kiefer and Vogelsang (2005).
The distribution of the test statistic (8) will depend on kernel (triangular in our case) and
the bandwidth chosen, which is set by default equal to T 0.5 , as suggested by Coroneo and
Iacone (2015). The results can be very different than those resulting from the traditional
asymptotic theory, where the test statistic would have the same distribution under the
null independently of the kernel and the bandwidth used.
Encompassing Test
Independently of whether the null hypothesis E[e2t − ĕ2t ] = 0 is rejected or not, it is
relevant to understand to what extent our model encompasses all the relevant information
of the benchmark, and the other way around. Because of the obvious symmetry of both
statements, we consider only the first one. If our forecast yt|Fi encompasses a given
benchmark y̆t|Fi , the difference between those benchmark forecasts and ours will not be
a relevant factor in explaining our own forecast error. In other words, the regression
coefficient λ will not be significantly different from zero in the following regression:
yt − yt|Fi = λ (y̆t|Fi − yt|Fi ) +ξt (9)
| {z } | {z }
et et −ĕt
m
yt = λy̆t|Fi +(1 − λ)yt|Fi + ξt (10)
Following Harvey, Leybourne and Newbold (1997), the statistical significance of the λ
coefficient in regression 9 can be used to reject the null hypothesis that our model encom-
passes the benchmark. In this case of rejection, equation 10 suggests that a combination
of the two forecasts would yield a more informative forecast.
By construction, the value of the coefficient of a regression e˘t = α(e˘t − et ) + ξt is equal
to 1 − λ, but it is not necessarily true that the rejection of the null hypothesis in the first
case implies the acceptance of the symmetric statement.
The test-statistic is computed as follows. When the null hypothesis is that our model
encompasses the benchmark, we define the sequence {dt }Tt=1 , where dt = et (et − e˘t ), and
5
d¯
we compute E1 = s , exactly as in equation 8.
2π fˆd (0)
T
Efficiency: Bias Test
In order to assess whether our forecasts are unbiased, we will simply test the statistical
significance of the average error. In some cases, the time series of forecast errors {et }Tt=1
may be autocorrelated to some extent even when they are based on a model with IID
innovations. In such cases, the variance associated to the estimate of the average forecast
error may be large. The test statistic has exactly the same form as the previous tests
discussed so far.
Efficiency: Autocorrelation Test
We will test here a second necessary condition for our forecasts to be efficient: absence
of autocorrelation. In the same spirit as the tests described above, we will assess the
statistical significance of the forecast errors’ autocorrelation. Thus, our sequence {dt }Tt=1
will be defined with dt = et et−1 .
Implementation in J Demetra+
Structure of the library
The code is structured as follows:
1. The class AccuracyTests contains all functions required to compute the statistic
defined in equation (8). It also incorporates instructions regarding how to calculate
the pvalues when depending on whether one wants to use standard asymptotics or
fixed-smoothing asymptotics (the input AsymptoticsType is required).
2. The class AccuracyTests is extended by each one of the classes containing the
tests listed in Table (A.1). The constructor of each one of these classes can generate
the tests when either the forecasts or the forecast errors are given as an input. These
classes also define the method to calculate the loss function dt , which is specific to
each test.
6
3. The class GlobalForecastingEvaluation can also be used to generate all
the tests. It contains BiasTest, EfficiencyTest, DieboldMarianoTest
and EncompassingTest and the corresponding methods to get the test statistics
and p-values described in Table A.1.
4. The class ForecastEvaluation contains methods to quantify errors: Root Mean
Squared Errors (RMSE), relative RMSE, Mean Absolute Errors (MAE), etc...
A simple example
Suppose we want evaluate the forecast of a model (ftm ) and compare them with those of
a benchmark (ftb ). The following points explain all the steps followed in the code below
to run all the tests:
• First we need to initialize the two competing forecast (i.e. benchmark vs model), all
the statistics we are going to calculate (RMSE, bias, autocorrelation, and encom-
passing weights) and the p-values corresponding to each one of the tests.
• Second, we initialize the eval object of the class GlobalForecastingEvaluation,
which will contain all test results. The inputs needed to run the tests are three time
series (our model’s forecasts, those of the benchmark, and the actual data, which is
the target) and the kind of distribution of the various test statistics under the null,
which is given by a normal distribution when
[Link].STANDARD_FIXED_B
is used.
• By choosing the option
[Link].HAR_FIXED_B
the distribution tabulated by Kiefer and Vogelsang (2005) is used.
• For each type of test, the bandwidth used to estimate the variance needs to be
specified. Otherwise, the default value will be used (T 1/2 ). The relevant statistics
for each test as well as the pvalues are obtained with a simple get command. Notice
7
that getPValue(twoSided) uses the logical argument true in order to get the
p-values of the two-sided test.
public void example() {
TsData[] series = {benchmark, model, target};
boolean twoSided = true;
double rmse = new double ;
double dmPval = new double ;
double bias = new double ;
double biasPval = new double ;
double arcorr = new double ;
double arPval = new double ;
double m_enc_bench = new double ;
double m_enc_bench_Pval = new double ;
double bench_enc_m = new double ;
double bench_enc_m_Pval = new double ;
// squared root of T
int bandwith = (int) [Link]([Link](), 1.0 / 2.0);
GlobalForecastingEvaluation eval = new GlobalForecastingEvaluation(model, benchmark,
target,
[Link].HAR_FIXED_B);
[Link]().setBandwith(bandwith);
dmPval = [Link]().getPValue(twoSided);
ForecastEvaluation feval = new ForecastEvaluation(model, benchmark, target);
rmse = [Link]();
[Link]().setBandwith(bandwith);
bias = [Link]().getAverageLoss();
biasPval = [Link]().getPValue(twoSided);
[Link]().setBandwith(bandwith);
arcorr = [Link]().calcCorelation();
arPval = [Link]().getPValue(twoSided);
[Link]().setBandwith(bandwith);
m_enc_bench = [Link]().calcWeights();
8
m_enc_bench_Pval = [Link]().getPValue(twoSided);
bench_enc_m = [Link]().calcWeights();
bench_enc_m_Pval = [Link]().getPValue(twoSided);
}
Example in the context of nowcasting
Let’s use a more complex example taken from Basselier, de Antonio and Langenus (2017)
in the context of a model for nowcasting, where blocks of data releases are used to update
the predictions. The aim is to determine which blocks of releases lead to significant
improvements of the forecasting accuracy.
The code has been reproduced below. The first variable series is an array of time
series, which includes the forecasts of an ARIMA model followed by the nowcasts from a
dynamic factor model obtained at different points in time, ranging from 90 days before the
end of the quarter (DFMPRE 90) to 44 days after the end of the quarter (DFMPOST 44).
The variable FLASH corresponds to the GDP growth, which is the target, or the variable
the model aims to predict.
The following code uses a simple loop to produce the test results for the eleven DFM
updates. For each one, we not only test the bias and autocorrelation, but we also check
whether it adds any value with respect to the previous forecast by means of both the
Diebold-Mariano and Encompassing tests. For example, is DFMPRE 90 better than
ARIMA? Is DFMPRE 75 better than DFMPRE 90? More precisely, as explained above,
the Diebold-Mariano test simply evaluates whether the difference in accuracy for each
pair is statistically significant, while the encompassing test checks whether the updated
forecast can be explained by its deviation with respect to the old one. Rejecting this
hypothesis implies that the new forecast encompasses the old one (i.e. the weight λ in
equation 10 above is not significantly different from zero).
public void testingUpdates() {
TsData[] series = {ARIMA, DFMPRE90, DFMPRE75, DFMPRE60, DFMPRE45, DFMPRE30, DFMPRE15,
DFM0, DFMPOST15, DFMPOST30, DFMPOST42, DFMPOST44, FLASH};
boolean twoSided = true;
9
double[] rmse = new double[11];
double[] dmPval = new double[11];
double[] bias = new double[11];
double[] biasPval = new double[11];
double[] arcorr = new double[11];
double[] arPval = new double[11];
double[] m_enc_bloom = new double[11];
double[] m_enc_bloom_Pval = new double[11];
double[] bloom_enc_m = new double[11];
double[] bloom_enc_m_Pval = new double[11];
for (int i = 1; i < [Link]; i++) {
// squared root of T
int bandwith = (int) [Link](series[i].getObsCount(), 1.0 / 2.0);
GlobalForecastingEvaluation eval = new GlobalForecastingEvaluation(series[i], series[i -
1], FLASH, [Link].HAR_FIXED_B);
[Link]().setBandwith(bandwith);
dmPval[i] = [Link]().getPValue(twoSided);
ForecastEvaluation feval = new ForecastEvaluation(series[i], ARIMA, FLASH);
rmse[i] = [Link]();
[Link]().setBandwith(bandwith);
bias[i] = [Link]().getAverageLoss();
biasPval[i] = [Link]().getPValue(twoSided);
[Link]().setBandwith(bandwith);
arcorr[i] = [Link]().calcCorelation();
arPval[i] = [Link]().getPValue(twoSided);
[Link]().setBandwith(bandwith);
m_enc_bloom[i] = [Link]().calcWeights();
m_enc_bloom_Pval[i] = [Link]().getPValue(twoSided);
bloom_enc_m[i] = [Link]().calcWeights();
bloom_enc_m_Pval[i] = [Link]().getPValue(twoSided);
[Link]("RMSE" + "\t" + rmse[i] + "\t" + dmPval[i] + "\t" + "Bias" + "\t" +
bias[i] + "\t" + biasPval[i] + "\t" + "CORR" + "\t" + arcorr[i] + "\t" + arPval[i] +
"\t" + "Weight on Update" + "\t" + (1 - m_enc_bloom[i]) + "\t" + "M_enc_bench Pval" +
"\t" + "Bench_enc_M Pval" + "\t" + bloom_enc_m_Pval[i]);
}
}
10
The outcome of theses tests has been organized in Table A.2. It includes bias, au-
tocorrelation, RMSE and the λ coefficient defined above, which is the weight given to
a benchmark forecasts that competes with their model’s. Statistical significance is high-
lighted with shades. Grey shaded areas in column FS-DM demonstrate which news blocks
have induced a significant change in the RMSE of the model, i.e. the null of equal accu-
racy between old (O) and updated (U) forecasts is rejected. The outcome of the DM test
may be considered jointly with the results of the encompassing tests. For a certain news
block to be considered relevant, the corresponding nowcasting update (U) should hold a
larger amount of information than the older nowcast (O) based on the previous informa-
tion set, while the old nowcast does not incorporate any useful information absent in the
new update. The last two colums of the table show that this is generally the case, with
some exceptions. That is, the null U encompases O is not rejected while O encompasses
U is rejected.
11
Table A.2: Statistical significance of each update based on fixed-smoothing (FS) asymp-
totics
Evaluation period: 2007.Q1 - 2015.Q1, T=25
FS-Efficiency FS-DM FS-Encompassing
(U)pdate vs (O)ld
Real-Time Updates bias corr RMSE U enc O O enc U
ARIMA -0.27 0.50 - - -
DFM -90 (d)ays -0.22 0.41 0.68 0.60 0.39
DFM -75 d -0.19 0.47 0.55 -0.60 1.59
DFM -60 d -0.12 0.55 0.52 0.26 0.51
DFM -45 d -0.14 0.54 0.54 1.48 -0.54
DFM -30 d -0.08 0.58 0.50 -0.20 1.07
DFM -15 d -0.13 0.46 0.41 -0.65 1.59
DFM 0 d (end of quarter) -0.06 0.45 0.38 -0.13 0.82
DFM +15 d -0.09 -0.11 0.27 -0.02 1.01
DFM +30 d -0.07 -0.08 0.26 -0.39 1.23
DFM +42 d -0.10 -0.06 0.26 0.27 0.66
DFM +44 d -0.06 -0.18 0.23 -0.17 1.03
Note: The FS-Efficiency multicolumn of his table reports bias and autocorrelation for the
forecast errors obtained at different horizons. The FS-DM and FS-Encompassing blocks
should be considered simultaneously. They aim to determine for each forecasting update (U)
whether there is any added value with respect to the old/last available forecast (O). The null
hypothesis of the Diebold-Mariano (DM) test is rejected when the difference in the squared
errors of U and O is significantly different from zero. For the two encompassing tests, the null
hypothesis states that the updated forecast (U) encompasses all the relevant information from
the old forecast (O) (or vice versa). When the null hypothesis can be rejected, this implies
that U can be improved by combining it with O. The combination weight associated to O (or
U ) is therefore reported below the “U enc O” test. In order to assess the added value of the
updated forecast, the DM null of equal forecast accuracy should be rejected and at the same
time the null “U enc O” and “O enc U” should be, respectively, not rejected and rejected.
Given the small size of our evaluation sample and the time-series correlation patterns, we
determine significance at the 5% , 10% and 20% level using the fixed-smoothing (FS)
asymptotics, as proposed by Coroneo and Iacone (2015).
12
References
[1] Coroneo, L. and F. Iacone (2016). “Comparing predictive accuracy in small samples
using fixed-smoothing asymptotics”. Discussion Papers University of York, 15/15.
[2] Basselier R., D. de Antonio Liedo, and G. Langenus (2015). “Nowcasting Euro Area
GDP: Understanding the Role of Qualitative Surveys”. , , x-X
[3] Diebold, F.X. and R.S. Mariano (1995). “Comparing Predictive Accuracy”. Journal
of Business and Economic Statistics, 13, 253-265.
[4] Kiefer, N.M., and T.J. Vogelsang (2005). “ A new asymptotic theory for
heteroskedasticity-autocorrelation robust tests”. Econometric Theory, 21, 1130-1164.
[5] Newey, W.K. and K.D. West (1987). “A simple, positive semi-definite, heteroskedas-
ticity and autocorrelation consistent covariance matrix”. Econometrica, 55(3), 703-
708.
13