Hypothesis: time series has an inverted-U shape.

How do we test this numerically?

My idea is to take the first difference of the variable and fit a linear model using the differentiated variable as endogenous variable and the time variable as exogenous.

$\Delta y_t = \beta_1 + \beta_2 t_t + \epsilon_t$

If the hypothesis is true, $\beta_2$ should be significantly less than zero.

If we try this approach with computer generated data, it can be seen that it works well.

enter image description here

lm(formula = dy ~ dt)

       Min         1Q     Median         3Q        Max 
-1.219e-15 -2.520e-16 -2.218e-17  1.827e-16  1.241e-15 

              Estimate Std. Error    t value Pr(>|t|)    
(Intercept)  1.210e-01  1.118e-16  1.082e+15   <2e-16 ***
dt          -2.000e-03  2.245e-18 -8.910e+14   <2e-16 ***
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 4.988e-16 on 82 degrees of freedom
Multiple R-squared:     1,  Adjusted R-squared:     1 
F-statistic: 7.939e+29 on 1 and 82 DF,  p-value: < 2.2e-16 

However, if a slight amount of noise is added to the data, this method falls apart catastrophically:

enter image description here

lm(formula = dy ~ dt)

     Min       1Q   Median       3Q      Max 
-0.96480 -0.21802  0.00826  0.24701  0.93200 

             Estimate Std. Error t value Pr(>|t|)
(Intercept)  0.114305   0.087548   1.306    0.195
dt          -0.001922   0.001758  -1.093    0.277

Residual standard error: 0.3907 on 82 degrees of freedom
Multiple R-squared: 0.01437,    Adjusted R-squared: 0.002345 
F-statistic: 1.195 on 1 and 82 DF,  p-value: 0.2775 

So, what is the alternative?


R code to generate the series and the plots:

t <- 1:85
y <- 0.12 * t - 0.001 * t^2 + rnorm(length(t), sd=0.25)
dt <- tail(t, -1)
dy <- tail(y, -1) - head(y, -1)

plot(t, y, ylim=c(-0.5, 4), pch=19, col='navy')
points(dt, dy, pch=19, col='purple')
legend(x=3, y=3.5, c('y','first difference'), pch=19, col=c('navy','purple'))

summary(lm(dy ~ dt))
  • $\begingroup$ Why don't you post the series with the noise and I will respond. $\endgroup$ – IrishStat Dec 13 '12 at 19:20
  • $\begingroup$ @IrishStat Done. $\endgroup$ – Ernest A Dec 13 '12 at 19:41
  • 1
    $\begingroup$ You could directly test the significance of a squared term with summary(lm(y~t+I(t^2)) (taking a one-sided test, since you apparently know beforehand that the series is concave). $\endgroup$ – Stephan Kolassa Dec 13 '12 at 21:20
  • 1
    $\begingroup$ @StephanKolassa The model that you suggest has colinearity problems, no? $\endgroup$ – Ernest A Dec 14 '12 at 15:28
  • 2
    $\begingroup$ @Ernest That's an interesting point. Using orthogonal polynomials cures the problem. E.g., u <- poly(t,2); summary(fit <- lm(y~u[,1]+u[,2])); library(car); vif(fit) shows the variance inflation factors are as low as possible (both equal to $1$). $\endgroup$ – whuber Dec 14 '12 at 20:50

If you use lm then you should check the residuals to see if they are autocorrelated or not. I guess they are not uncorrelated and hence your t-test are not valid (this is true also for the case of summary(lm(y~t+I(t^2)). This is basiacally beacuse there is a time variable involved in your lm.

I recommend to use Generalized Least Square approach in order to test the quadratic effect and take into account the autocorrelated problem. For example if you assume the autoregressive of order two (see below) for the residuals of your lm (i.e. $e_t=\phi_1 e_{t-1}+\phi_2 e_{t-2}+\nu_t$, where $\nu_t$ is white noise), then the code would be like


Note: You should model the error terms correctly first (i.e. finding the order of $p$ and $q$) maybe by ckecking the ACF or PACF of the residuals in your lm. In above, I assumed AR(2). More complicated ARMA model can be considered and tested.

  • $\begingroup$ +1: Thank you for illustrating the application of GLS. In your example, to match the code in the question, the argument data=your.data., should be removed. It's also a good idea to use orthogonal polynomials rather than direct computation of the square of t. $\endgroup$ – whuber Dec 14 '12 at 20:54
  • $\begingroup$ YW, I removed the argument data in the code. $\endgroup$ – Stat Dec 14 '12 at 21:26

Your data:

enter image description here

If you know a priori what the form of the equation is then as others have pointed out it is trivial to estimate parameters. A more general solution is to characterize the data with an ARIMA model that may or may not include pulses, level shifts, local time trends, changes in parameters over time and/or changed in error variance over time. This approach enter image description here

yields an equation of the form

enter image description here

The residuals suggest an adequate model

enter image description here

Remember all models are wrong but some models are useful (G.E.P.Box). The final model summary reports a mean-square error of the errors to be .1 which should agree with your simulation

enter image description here

In terms of transparency, I am one of the developers of AUTOBOX, which I used here to illustrate what could be done to form a possible model/prediction.

  • 2
    $\begingroup$ Re your comments about wrongness and usefulness of models: It is striking that you report an eight-parameter model for a (mathematically and statistically simple) process that was generated with only three parameters. It would be natural to characterize that as overfitting. How would you respond? $\endgroup$ – whuber Dec 14 '12 at 19:28
  • $\begingroup$ @huber one could always restrict the number of pulses that were identifiable or set a more stringent alpha level for retention. You state that only three parameters were used to form the signal but the error process introduced some possible irregularities.The 5 pulses were attributed to the irregularities in the added noise (potential over-parameterization without penalty) while the double difference and the AR(2) constituted the equivalence of 4 coefficients in the signal/prediction equation. What approach would you suggest to separate the data to signal and noise ? $\endgroup$ – IrishStat Dec 14 '12 at 20:18
  • 1
    $\begingroup$ When you look at the code that generated these data, you see that the "error process" is iid normal. (One of the three parameters I referred to is its variance.) So it seems you're saying that your method has a tendency to identify many "pulses" and "level shifts" where in fact those do not really happen at all (certainly not at the rate of 5 in 81 observations). A more nuanced understanding of the "wrong" versus "useful" quotation will bring in the idea parsimony, which seems to be lacking here. $\endgroup$ – whuber Dec 14 '12 at 20:26
  • $\begingroup$ The error process may be iid normal but the defaault is set at alpha of .05 which may may incorrctly identify some anomnalies in the data. To remedy this one cam always set a maximumum for the number of identified pulses, seasonal pulses, level shifts , local time trends OR set an alpha level smaller say .01 to restrict the "exhuberance". It is all controllable by you as the software allows this kind of user-specification if you so desire and it appears you desire to do so for purposes of parsimony. $\endgroup$ – IrishStat Dec 14 '12 at 20:38
  • $\begingroup$ Thank you for clarifying that--it's good to know that the analysis can be controlled in this way. $\endgroup$ – whuber Dec 14 '12 at 20:55

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.