Fitting NLS algorithm

Question

I have the list of time series. I am fitting these series with the formula $$y=ax^2\exp(b*x)$$. It should be noted that parameter $b$ in the formula must be negative as this reflects the behaviour of time series. An example plot for the fitted regression output using nls() in R is:

Initially, I was faced with the problem of choosing start values for $a$ and $b$. I chose such values so that the fitted value evaluated at the maximum for $y$ is the same.

However with some time series I faced with the next problem. Let's consider this time series:

When applying my methodology, nls() gave me a positive value of $b$ which is unsatisfactory for me although the fitted values were highly concordant. However when I tried to introduce boundaries in the nls() function such as upper(a=10,b=-0.0001), I found that no matter what value of $b$ I put, the resulting value of $b$ does not converge to any stable value and the fitted curve is highly discordant. Numerically, the parameters were consistent with what I expected to find, but the fitted values were uncalibrated.

Can anybody give any suggestions of how can I deal with such time series without changing the general formula? Thanks!

Answer 1

This is clearly an issue of an optimizer yielding values on the boundary of your parameter space. I suspect if you obtained the Hessian from nls(), it would be singular. One can see straight away that the curve you displayed in the second figure is best fit by $b = 0$ exactly and $a$ scaling the quadratic curve. In general, it is a very bad idea to enforce constraints on parameters with general solvers such as this one. You've already suggested that the wavelet functional form thus specified is immutable, but if it were me running the analysis, I would strongly consider an alternate parameterization:

$$ y = ax^2 \exp \left( -\left( \frac{x}{b} \right) ^2 \right) $$

Note additionally the issue involved with specifying starting values. As I mentioned, the second graphic strongly suggests a fully quadratic trend. This is especially driven by your method of choosing starting values. I've not even cranked the numbers, but I would suspect you solved for $a \approx 1/250^2$ and $b \approx 0$. For very large $x$, the exponential term dominates the quadratic term, hence $b$ is forced toward 0 to accommodate trends that show large non-zero values in the tails. So my first recommendation is to use another rule for choosing starting values. I don't see why $a^{(0)},b^{(0)}=1$ is out of consideration.

Consider additionally a change of variable. $Y$ is positive valued, I assume. That means you've violating some grave assumptions we would have to exploit to fully utilize nls. Namely, the assumption that errors are symmetric around the functional mean. If you simply log transform $y$, you'll have much better convergence values, and it's possible to transform the coefficients to their original scale ($b$ is consistently estimated, whereas estimates of $log(a)$ must be exponentiated).

$$ \log(y) = 2\left(\log(a) + \log(x) \right) + bx$$

This is assuming of course that there are errors in the data you've measured. Otherwise, some type of wavelet decomposition would be better suited to modeling what appear to be signals.