I'm having trouble understanding the way that Gretl computes forecasts for ARIMA(0,0,1)(0,0,1) models. If I generate a random data set and fit this model (using either algorithm) to get the estimates theta_1 and Theta_1, and then generate a forecast, then I would expect the forecast to be computed using the polynomial 1 + theta_1*L + Theta_1*L^4 + theta_1*Theta_1*L^5 [... But] after a little experimenting, I have found that if I use the polynomial 1 + theta_1*L + Theta_1*L^4 - theta_1*Theta_1*L^5 I get exactly the same forecast as Gretl.