How to detect specific behavior in time series?

Question

I was not quite sure what the right SE for this was, so I posted this also here on SciComp. 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:

1. a power-law (the $(t_f-t)^\gamma$-term)
2. 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:

1. smooth the signal to get rid of the noise and the periodic component
2. take the first two derivatives numerically.
3. Periods of consecutive positive first derivative are a necessary condition for the given power law -- continue analyzing only these periods ($t_0$ estimated)
4. The division of the first derivative by the second makes it possible to estimate $t_f$.
5. 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.
6. 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:

1. Take an initial $t_f^{\text{trial}} > t_B$, so some time after the beginning of the time series.
2. Compute $\ln(t_f^{\text{trial}}-t)$, from $t=t_B$ to just before $t_f^{\text{trial}}$.
3. 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).
4. Determine all the peaks in the frequency domain and save them.
5. Repeat from the top with $t_f^{\text{trial}} \leftarrow t_f^{\text{trial}} + \Delta t$
6. 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.
7. For all peaks and $t_f$ thus found, find the corresponding interval that maximizes the power in the peak ($t_0$ estimated).
8. 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:

1. The estimates for $t_0$ and $t_f$ from both approaches will generally differ. determine the smallest overlapping interval $(t_0, t_f)$.
2. Throw the data from this interval, as well as all the initial estimates, into a non-linear least squares fitter to improve the fit.
3. 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:

1. seems to be extremely sensitive to noise
2. I don't know how to determine automatically what a "good enough" smoothing is
3. 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:

1. 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
2. detecting which peak correspond to which other peak from one estimate for $t_f$ to the next, is rather difficult to automate.
3. 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?

Answer 1

An approach that might be worthwhile investigating is some sort of on-line method. Because of the highly non-linear nature of your signal model, the standard Kalman filter is out, but an Extended Kalman filter (or, perhaps, a particle filter) might be worth looking at.

The idea is that you set up your signal model state to be: $$ x_t = \left[ \alpha\ \beta\ \gamma\ \delta\ \phi\ \omega\ t_f \right]^T $$ with appropriate possible changes to the parameters modeled by the process noise covariance: $$ x_{t+1} = x_t + w_k $$ with $Q_k$ being the covariance of $w_k$.

The output equation would be where all the nonlinearities are encoded.

The trick will be to see how the (Extended) Kalman filter's error covariance differs between when you have this signal present, and when you don't. The CUSUM algorithm applied to some statistic you will need to define might work. See Basseville and Nikiforov for details of that and other change detection algorithms.

It's not exactly the same thing, but this paper suggests that the approach may be feasible. Their first example has a somewhat nonlinear system: