I was not quite sure what the right SE for this was, so I posted this also here on DSP. Please tell me which one to remove :)
Problem statement
I have a few hundred unrelated time series, say $P_i(t)$, with $i=1,2,...,N$. The sample times are all unequally spaced.
It is known that during certain (unknown) intervals, the trend in some of the time series is well-approximated by
$$ P_i(t) \approx \alpha - \beta \left(t_f-t\right)^\gamma\left(1 + \delta\cos\left(\omega\ln\left(t_f-t\right)+\phi\right)\right) $$
with
$$ t_0 \leq t < t_f\\ \alpha,\beta > 0\\ 0\leq\gamma,\delta\leq1\\ 0\leq\phi<2\pi $$
and $\omega>0$ usually "small", but not really restricted.
An example of what such a function looks like:
After this interval ($t > t_f$), there is nearly always a rapid decline, after which the signal continues as it did before the interval ($t<t_0$). Of course, each series also has zero-mean, normally distributed noise of varying strength superimposed on it, hence the approximate sign.
My task:
- to determine the intervals where this behavior occurs,
- find best-estimates for all the parameters listed above.
My first stab at it
Since each of these time series can be quite long and voluminous, and I have quite many series to analyze (and that number will likely grow over time), visual inspection is out of the question. Moreover: often, it is quite easy to miss this pattern by just looking at the signal.
The pattern consists of two main components:
- a power-law (the $(t_f-t)^\gamma$-term)
- an oscillation with logarithmically varying period (the $\cos\left(\omega\ln\left(t_f-t\right)\right)$-term)
To detect the power law, I thought of the following procedure:
- smooth the signal to get rid of the noise and the periodic component
- take the first two derivatives numerically.
- Periods of consecutive positive first derivative are a necessary condition for the given power law -- continue analyzing only these periods ($t_0$ estimated)
- The division of the first derivative by the second makes it possible to estimate $t_f$.
- The standard deviation and mean of each element's estimate for $t_f$, as well as the demand that $t_f > t$ and $t < t_E$ (the final time at which the first derivative is positive), give a measure for the reliability for this estimate.
- If it is reliable enough, it is fairly straightforward to backtrack everything and come up with estimates for $\gamma$, $\beta$ and $\alpha$ (in that order).
To detect the oscillating component, I came up with the following:
- Take an initial $t_f^{\text{trial}} > t_B$, so some time after the beginning of the time series.
- Compute $\ln(t_f^{\text{trial}}-t)$, from $t=t_B$ to just before $t_f^{\text{trial}}$.
- Compute the Fourier transform of the transformed time series $P_T(t) = P_i(\ln(t_f^{\text{trial}}-t))$ (via Lomb/Scargle, because the new sample times are unequally and logarithmically spaced).
- Determine all the peaks in the frequency domain and save them.
- Repeat from the top with $t_f^{\text{trial}}\leftarrow t_f^{\text{trial}} + \Delta t$
- The progression of the peaks will roughly follow a parabola-shaped path in the frequency domain if there is a periodic component somewhere. The "right" $t_f$ will be found when the maximum power in the frequency has been found.
- For all peaks and $t_f$ thus found, find the corresponding interval that maximizes the power in the peak ($t_0$ estimated).
- with all this information, it is fairly straightforward to come up with initial estimates of $t_f$, $\omega$ and $\phi$.
Then, these two approaches are to be combined:
- The estimates for $t_0$ and $t_f$ from both approaches will generally differ. determine the smallest overlapping interval $(t_0, t_f)$.
- Throw the data from this interval, as well as all the initial estimates, into a non-linear least squares fitter to improve the fit.
- Compute a couple of measures related to the goodness of fit, which will accompany the results.
Why I'm looking for another approach
Well, that all sure sounds very nice, but of course I wouldn't be asking a question here if it all worked as well as it sounds :)
The method I use to detect the presence and parameters of the power law:
- seems to be extremely sensitive to noise
- I don't know how to determine automatically what a "good enough" smoothing is
- All smoothing algorithms I have tried are not very good at removing the periodic component, throwing all the estimates way off. Moreover, the (log-)periodic component can make the derivatives negative.
I could take the log of the data and detect linearity, but that suffers from the same problem -- the periodic component (as well as the unknown offset $t_f$) seems to make that unreliable.
The method I use to detect the periodic component:
- is very computationally intensive; Lomb/Scargle is certainly not as fast as an FFT would be. And as I mentioned earlier, each time series can be quite long, so the number of times the Fourier transform needs to be computed can also be quite large
- detecting which peak correspond to which other peak from one estimate for $t_f$ to the next, is rather difficult to automate.
- transforming the sample times into (roughly) logarithmically-spaced sampling times makes it very hard to detect the start/end of the interval with a decent accuracy.
I'm kind of stuck, and I need some new inspiration. Any suggestions?