# Fitting NLS algorithm

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!

• When the best fit has $b\gt 0$, if you want to constrain $b\le 0$ you can set its value manually. For instance, if you set $b=0$ you will get a quadratic fit. For smaller values of $b$ you must set them manually (according to some arbitrary criterion). Incidentally, to find decent starting values simply regress $\log(y)-2\log(x)$ against $x$: the intercept estimates $a$ and the slope estimates $b$.
– whuber
Commented Nov 20, 2013 at 17:31
• aslo consider trying - minpack.lm -. I'm not sure that it is justified, but I find the optimizer converges much more reliably than the base nls Commented Nov 20, 2013 at 18:48

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.

• +1 This is a really good answer. I'm sorry it hasn't been noticed before.
– whuber
Commented Apr 28, 2014 at 21:56