The commonly used RMSE (Fox, 1981) quantifies the differences between predicted and observed values, and thus indicates how far the forecasts are from actual data. A few major outliers in the series can skew the RMSE statistic substantially because the effect of each deviation on the RMSE is proportional to the size of the squared error. The overall, nondimensional measure of the accuracy of forecasts MASE (Hyndmann and Koehler, 2006) is less sensitive to outliers than the RMSE. The MASE is recommended for determining comparative accuracy of forecasts (Franses, 2016) because it examines the performance of forecasts, relative to a benchmark forecast. It is calculated as the average of the absolute value of the difference between the forecast and the actual value divided by the scale determined by using a random walk model (naïve reference model on the history prior to the period of data held back from model training). MASE<1 indicates that the forecast model is superior to a random walk.

The correlation coefficient between estimates and observations (Addiscott and Whitmore, 1987), -1 (anti-correlation) -1≤R≤1 (perfect correlation), assesses linear relationships, in that forecasted values may show a continuous increase or decrease as actual values increase or decrease. Its extent is not consistently related to the accuracy of the estimates.

Forecasts online means (http://smoothforecast.com/SmoothForecast/index.jsp) were used to employ simulation model with support of Excel spreadhseet. The statistics were assessed interactively using Statistics Software STATGRAPHIC Online (http://www.statgraphicsonline.com) and WESSA R–JAVA web (Wessa, 2012).

In order to quantify long-range dependence and appraise the cyclical-trend patterns in the series, we estimated the Hurst (H ) exponent (rate of chaos), which is related to the fractal dimension (D =2-H ) of the series. Long memory occurs when 0.5<H <1.0, that is, events that are far apart are correlated because correlations tend to decay very slowly. On the contrary, short-range dependence (0.0<H <0.5) is characterized by quickly decaying correlations, i.e. past trends tend to revert in the future (an up value is more likely followed by a down value). Calculating the Hurst exponent is not straightforward because it can only be estimated and several methods are available to estimate it, which often produce conflicting estimates (Karagiannis et al., 2002). Using SELFIS (http://www.cs.ucr.edu/~tkarag/Selfis/Selfis.html), we have referred to two methods, which are both credited to be good enough to estimate H (Belov et al., 2006): the widely used (Yin et al., 2009) rescaled range analysis (R/S method) and the ratio variance of residuals method, which is known to be unbiased almost through all Hurst range (Sheng and Chen, 2009). Long-memory in the occurrence of PDSI values was also analysed to see if the memory characteristic is correlated with the length of the time series. To determine whether this characteristic changes over time, the Hurst exponent was not only estimated for the full time series (1801-1914), but also for a shorter series starting in 1901 (the most recent period, which is also the period held out of the calibration process).

An alternative methodology based on transfer function models was also implemented to compare two methodologies: seasonal exponential smoothing that uses the temporal dependence structure of the time series itself to reproduce the time series behavior in the future, versus the linear transfer function approach where input time series potentially impacting the drought behavior at large spatial scales are used as explanatory time serie variables in a lagged regression model. The methodology called transfer fuction models, was introduced by Box and Jenkins (1970) and re-visited by Shumway and Stoffer (2016) to model the impact of El Niño on fish recruitment.

In a transfer function model the output series Y(t) (in this case PDSI) can be represented as:

Y(t)=α1(B) X1(t)+α2(B)X2(t)+…+αk(B)Xk(t) + η(t)

where X1(t), X2(t), …, Xk(t) are the input time series to be considered as explanatory variables contributing to the temporal dynamics of the output series Y(t). The terms α1(B), α2(B), …, αk(B) are fractional polynomials in the backward-shit operator B (such that Bs(X(t))= X(t-s)) of the form:

αi(B)=δi(B) Bdi/ωi(B)

where δi(B)= (δ0i+δ1i B1 +…+δsiBsi ) and ωi(B)= (1-ω1i B1-…-ωri Bri) and si, ri and di are the orders of the transfer function for each input time series Xi(t). The numerator polynomial (order si) controls the lagged effect on the input variable Xi(t) and the denominator (order ri ) controls the lagged effect of variable Y(t). The back-shift operator order di and the orders si and ri are estimated by calculating the cross-correlation function between the pre-whitened input series Xi(t) (using an ARIMA model, such that the filter Φ(B)/Θ(B) applied to Xi(t) produces a white noise ν(t)) and the filtered output series Y(t) using the same ARIMA filter as for the input series. Once the orders si, ri and di are identified, the corresponding linear transfer function is fitted using the R software. In our case i=1,2, since two imput series are used as covariates.