ETS models can produce bad roots of characteristic equation that cannot be refit

[mitchelloharawild]

Ref: https://stackoverflow.com/questions/68673653/why-am-i-getting-parameters-out-of-range-error-after-fitting-a-default-ets-mod

cc @robjhyndman - I think this is due a difference in checks of characteristic equation roots between the C++ during optimisation (used in model()/estimate()), and the R check during initialisation (used in refit())

fable/R/etsmodel.R

Lines 503 to 512 in 0f3c42f

	# End of easy tests. Now use characteristic equation
	P <- c(phi * (1 - alpha - gamma), alpha + beta - alpha * phi + gamma - 1, rep(alpha + beta - alpha * phi, m - 2), (alpha + beta - phi), 1)
	roots <- polyroot(P)
	# cat("maxpolyroots: ", max(abs(roots)), "\n")
	if (max(abs(roots)) > 1 + 1e-10) {
	return(0)
	}
	}

From R, we have P[2] = alpha + beta - alpha * phi + gamma - 1

fable/src/etsTargetFunction.cpp

Lines 280 to 306 in 0f3c42f

	// End of easy tests. Now use characteristic equation
	std::vector<double> op_new(m+2, alpha+beta-alpha*phi);
	// P <- c(phi*(1-alpha-gamma),alpha+beta-alpha*phi+gamma-1,rep(alpha+beta-alpha*phi,m-2),(alpha+beta-phi),1)
	op_new[0] = phi*(1-alpha-gamma);
	op_new[1] = phi*(alpha+beta-alpha*phi+gamma-1);
	op_new[m] = alpha+beta-phi;
	op_new[m+1] = 1;
	Environment base("package:base");
	Function polyroot = base["polyroot"];
	Function abs = base["abs"];
	NumericVector res = abs(polyroot(op_new));
	// Rprintf("alpha = %f, beta = %f, gamma = %f, phi = %f, m = %i\n",
	// alpha, beta, gamma, phi, m);
	// Rprintf("C: c(");
	// for(int i=0;i<opr.size();i++) {
	// Rprintf("%f, ", opr[i]);
	// }
	// Rprintf(")\n");
	// Rprintf("C_new: c(");
	//Rprintf("maxpolyroot: %f\n", max);
	double max_root = max(res);
	if(max_root > 1+1e-10) return(false);

From C++ we have op_new[1] = phi*(alpha+beta-alpha*phi+gamma-1)

Which is correct?

[robjhyndman]

The R code is correct (eq 10.5, p156 of the Springer book https://paperpile.com/shared/FoKgkF). You can see from the comment in the C++ code on line 283 that this was what was intended but somehow the phi multiple has been included. I guess this happened when the C++ version was introduced in the forecast package (cc @cbergmeir). It will need to be fixed there as well.

[cbergmeir]

Yes, seems Rob is correct and I made a mistake there when translating from R to C++...and then nobody noticed for 8 years. Sorry about that.