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.