Examples are from:

• Shumway, R. H. Stoffer D. S. (2011) Time Series Analysis and Its Applications With R Examples. New York: Springer-Verlag.
• Legal PDF copy of this book can be obtained from: http://www.stat.pitt.edu/stoffer/tsa3/tsa3.pdf
• Also, there are some issues with time series tools in base R (stats), particularly in the arima() function. See the following for some examples: http://www.stat.pitt.edu/stoffer/tsa3/Rissues.htm

The mean function is defined as: \[ \mu_{t} = E(x_{t}) = \int_{-\infty}^{\infty}xf_{t}(x)\ dx \]

The autocovariance function is the second moment: \[ \gamma(s, t) = cov(x_{s}, x_{t}) = E[(x_{s} - \mu_{s})(x_{t} - \mu_{t})], \quad \forall s,t \]

Weakly stationary process is

• The mean function \(\mu_{t}\) is constant and does not depend on time \(t\) and
• The autocovariance function \(\gamma(s, t)\) depends on \(s\) and \(t\) only through \(|s - t|\)

Random Walk with a drift:

\[ x_{t} = \delta t + \sum_{i=1}^{t}w_{i}\] \[ \mu_{t} = E(x_{t}) = \delta t + \sum_{i=1}^{t}E(w_{i}) = \delta t\]

set.seed(1)
n <- 200; delta <- 0.3; t <- 1:n
w <- cumsum(rnorm(n))
x_t <- delta * t + w

• Since the mean function, \(E(x_{t}) = \mu_{t}\) of stationary time series is independent of time \(t\), \(\mu_{t} = \mu\).
• If let \(s = t + h\) where \(t\) is the lag: \[ \gamma(t + h, t) = cov(x_{t+h}, x_{t} = x_{0}) = \gamma(h, 0) \equiv \gamma(h) \]
• For the stationary time seties, the autocovariance function becomes: \[ \gamma(h) = cov(x_{t + h}, x_{t}) = E[(x_{t+h} - \mu)(x_{t} - \mu)] \]
• And the autocorrelation function becomes: \[ \rho (h) = \frac{\gamma(t + h, t)}{\sqrt{\gamma(t + h, t + h) \gamma(t, t)}} = \frac{\gamma(h)}{\gamma(0)} \]

\[ x_{t} = \phi x_{t-1} + w_{t} \]

ar1 <- arima.sim(list(order = c(1, 0, 0), ar = 0.7), n = n)

ar1_fit <- arima(ar1, order = c(1, 0, 0), method = "ML")
ar1_fit

## 
## Call:
## arima(x = ar1, order = c(1, 0, 0), method = "ML")
## 
## Coefficients:
##         ar1  intercept
##       0.641      -0.17
## s.e.  0.054       0.21
## 
## sigma^2 estimated as 1.18:  log likelihood = -301,  aic = 607

Another way of writing the AR1 model: \[ y_{t} \sim Normal(\alpha + \phi y_{t-1}, \sigma) \]

data {
  int<lower=0> N;
  vector[N] y;
}
parameters {
  real alpha;
  real phi;
  real<lower=0> sigma;
} model {
  sigma ~ normal(0, 1);
  alpha ~ normal(0, 1);
  phi ~ normal(0, 1);
  for (n in 2:N)
    y[n] ~ normal(alpha + phi * y[n-1], sigma);
}

data <- list(N = n, y = ar1)
ar1_stan_fit <- stan("models/ar1.stan", data = data)
phi <- extract(ar1_stan_fit, 'phi')$phi
quantile(phi, probs = c(.025, .25, .50, .75, .975))

## 2.5%  25%  50%  75%  98% 
## 0.54 0.61 0.64 0.67 0.75

sigma <- extract(ar1_stan_fit, 'sigma')$sigma
quantile(sigma, probs = c(.025, .25, .50, .75, .975))

## 2.5%  25%  50%  75%  98% 
##  1.0  1.1  1.1  1.1  1.2

• Red = True AR1 coefecient
• Blue = MLE AR1 from arima()

rw <- cumsum(rnorm(100, 0, 1))
plot.ts(rw)

data <- list(N = 100, y = rw)
ar1_stan_fit <- stan("models/ar1.stan", data = data)
phi <- extract(ar1_stan_fit, 'phi')$phi
mean(phi > 1)

## [1] 0.013

data {
  ...
}
parameters {
  ...
  // bad idea
  real<lower=-1, upper=1> phi;
  ...
} model {
  ...
}

data <- list(N = 100, y = rw)
ar1_stan_fit <- stan("models/ar1_bad.stan", data = data)
phi <- extract(ar1_stan_fit, 'phi')$phi
mean(phi > 1)

## [1] 0