3
$\begingroup$

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:

enter image description here

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?

$\endgroup$
1
$\begingroup$

If there's some model for the signal for $t < t_0$ and $t > t_f$, perhaps it would be possible to try to fit the full signal using a Metropolis-Hastings or Goodman-Weare algorithm. This could then also yield some information about the probability distributions for the parameters involved.

$\endgroup$
  • $\begingroup$ Generally, there isn't a model outside the interval, or if there is, it is not generally known. Would that matter much for either algorithm you propose? $\endgroup$ – Rody Oldenhuis Nov 3 '13 at 9:08
  • $\begingroup$ Not really sure since I've only done fits in case the model for an entire dataset was known. But my gut feeling is that you would need to be able to specify at least something about the data in the other parts of the signal. Otherwise the algorithm may ignore the actually interesting part of the signal because it thinks it is in the t > t_f range. Do you perhaps have a more complete data set online somewhere? May help so make things more clear... $\endgroup$ – brm Nov 3 '13 at 11:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.