0

I have a time series consisting of a mean, a cosine harmonic wave, and seasonal ARIMA errors, and the seasonality period of both is 7 days:

$$ \begin{align} x_t &= \mu + h_t + u_t \\ h_t &= Amp\cdot \cos(\omega_s t - \varphi) \; \;\;\text{with}\;\; \omega_s = \frac{2\pi}{7}\\ u_t &= SARIMA(1,0,1)(1,0,0)[7] \end{align} $$

I nead to estimate SARIMA parameters ($\phi$, $\theta$, $\Phi$) as well as the amplitude ($Amp$) and phase shift ($\varphi$, which indicates where the harmonic element peaks) of $h_t$. Theoretically, one can model $h_t$ with a Fourier transform of order 1:

$$ f_t = a \cdot \cos(\omega_s) + b \cdot \sin(\omega_s) $$

By including $\cos(\omega_s)$ and $\sin(\omega_s)$ as an exogeneous regressor in a SARIMA model, and estimate $a$ and $b$ as regression coefficients with their standard errors. However, I need point estimates (and their respective standard errors) of $Amp$ and $\varphi$ instead of $a$ and $b$.

What I had in mind was to write down $h_t$ with the parameters of $f_t$ (as shown here, but with a sine) via

$$ h_t = \begin{cases} \sqrt{a^2 + b^2} \cdot \cos(\omega_s t - \arctan\frac{a}{b}) & b > 0 \\ a \cdot \cos(\omega_s t - \frac{\pi}{2}) & b = 0 \\ -\sqrt{a^2 + b^2} \cdot \cos(\omega_s t - \arctan\frac{a}{b}) & b < 0 \end{cases} $$

Thus

$$ \begin{align} Amp = &\begin{cases} \sqrt{a^2 + b^2} & b > 0 \\ a, \; \; & b = 0 \\ -\sqrt{a^2 + b^2} & b < 0 \end{cases} \\ \\ \varphi = &\begin{cases} \arctan\frac{a}{b} & b \neq 0 \\ \frac{\pi}{2} & b = 0 \end{cases} \end{align} $$

And approximately calculate the standard errors of $Amp$ and $\varphi$ (as nonlinear combinations of $a$ and $b$) using, e.g., the delta method as explained in this and this answers.

However, since $Amp$ is piecewise—and hence non-differentiable—the differentiability assumption of the delta method is violated. Apart from the standard errors, intuitively, I cannot imagine what would happen to $Amp$ and $\varphi$ if the confidence intervals of $a$ and/or $b$ covers (or gets close) to zero.

Is there a way around it? Thanks in advance!

MWE and reproducible code:

library(forecast)
library(sarima)

set.seed(1)
N <- 1000 # time series length


Mean


mu <- 1


Harmonic component of the mean structure


amp <- 3
omega_s <- 2 * pi/7 # angular frequency for f = 1/7
phase_peak <- 3.5 * omega_s


h_t <- amp * cos(c(1:N) - phase_peak)


Seasonal arima noise


u_t <- sarima::sim_sarima(n = N,
               model = list(ar = 0.5,
                            ma = 0.2,
                            sar = 0.3,
                            nseasons = 7,
                            sigma2 = 16),
               n.start = 100
               )


x_t <- mu + h_t + u_t
x_t <- ts(x_t, frequency = 7)


Making sine-cosine pairs for external regressors


rh <- forecast::fourier(x_t, K = 1)


Estimating the model


arima(x_t,
      order = c(1,0,1),
      seasonal = c(1,0,0),
      xreg = rh)
#> 
#> Call:
#> arima(x = x_t, order = c(1, 0, 1), seasonal = c(1, 0, 0), xreg = rh)
#> 
#> Coefficients:
#>          ar1     ma1    sar1  intercept    S1-7     C1-7
#>       0.4134  0.2434  0.3147     0.6642  0.4324  -3.1029
#> s.e.  0.0484  0.0518  0.0301     0.4062  0.3910   0.3909
#> 
#> sigma^2 estimated as 17.38:  log likelihood = -2847.3,  aic = 5708.6

Created on 2023-06-18 with reprex v2.0.2

Which shows good estimates of all parameters, and the manual calculation of the amplitude ($Amp = \sqrt{0.432^2 + (-3.102)^2} = 3.133$) and phase shift ($\varphi = \arctan \frac{-3.103}{0.432} = -1.286$) seem "reasonable".

psyguy
  • 311

0 Answers0