# Time Series Analysis

A time series is a sequence of data points, each associated with a time. In our example, we will work with a time series of daily temperatures in the city of Melbourne, Australia over a period of a few years. Let $x$ be the vector of the time series, and $x_i$ denote the temperature in Melbourne on day $i$. Here is a picture of the time series:

```
using Plots, Convex, ECOS, DelimitedFiles
const MOI = Convex.MOI
aux(str) = joinpath(@__DIR__, "aux_files", str) # path to auxiliary files
temps = readdlm(aux("melbourne_temps.txt"), ',')
n = size(temps, 1)
plot(
1:n,
temps[1:n],
ylabel = "Temperature (°C)",
label = "data",
xlabel = "Time (days)",
xticks = 0:365:n,
)
```

We can quickly compute the mean of the time series to be $11.2$. If we were to always guess the mean as the temperature of Melbourne on a given day, the RMS error of our guesswork would be $4.1$. We'll try to lower this RMS error by coming up with better ways to model the temperature than guessing the mean.

A simple way to model this time series would be to find a smooth curve that approximates the yearly ups and downs. We can represent this model as a vector $s$ where $s_i$ denotes the temperature on the $i$-th day. To force this trend to repeat yearly, we simply want

\[s_i = s_{i + 365}\]

for each applicable $i$.

We also want our model to have two more properties:

- The first is that the temperature on each day in our model should be relatively close to the actual temperature of that day.
- The second is that our model needs to be smooth, so the change in temperature from day to day should be relatively small. The following objective would capture both properties:

\[\sum_{i = 1}^n (s_i - x_i)^2 + \lambda \sum_{i = 2}^n(s_i - s_{i - 1})^2\]

where $\lambda$ is the smoothing parameter. The larger $\lambda$ is, the smoother our model will be.

The following code uses Convex to find and plot the model:

```
yearly = Variable(n)
eq_constraints = [yearly[i] == yearly[i-365] for i in 365+1:n]
smoothing = 100
smooth_objective = sumsquares(yearly[1:n-1] - yearly[2:n])
problem = minimize(
sumsquares(temps - yearly) + smoothing * smooth_objective,
eq_constraints,
);
solve!(
problem,
MOI.OptimizerWithAttributes(ECOS.Optimizer, "maxit" => 200, "verbose" => 0),
)
residuals = temps - evaluate(yearly)
# Plot smooth fit
plot(1:n, temps[1:n], label = "data")
plot!(
1:n,
evaluate(yearly)[1:n],
linewidth = 2,
label = "smooth fit",
ylabel = "Temperature (°C)",
xticks = 0:365:n,
xlabel = "Time (days)",
)
```

We can also plot the residual temperatures, $r$, defined as $r = x - s$.

```
# Plot residuals for a few days
plot(1:100, residuals[1:100], ylabel = "Residuals", xlabel = "Time (days)")
```

`root_mean_square_error = sqrt(sum(x -> x^2, residuals) / length(residuals))`

`2.7488719721107198`

Our smooth model has a RMS error of $2.7$, a significant improvement from just guessing the mean, but we can do better.

We now make the hypothesis that the residual temperature on a given day is some linear combination of the previous $5$ days. Such a model is called autoregressive. We are essentially trying to fit the residuals as a function of other parts of the data itself. We want to find a vector of coefficients $a$ such that

\[\text{r}(i) \approx \sum_{j = 1}^5 a_j \text{r}(i - j)\]

This can be done by simply minimizing the following sum of squares objective

\[\sum_{i = 6}^n \left(\text{r}(i) - \sum_{j = 1}^5 a_j \text{r}(i - j)\right)^2\]

The following Convex code solves this problem and plots our autoregressive model against the actual residual temperatures:

```
# Generate the residuals matrix
ar_len = 5
residuals_mat = Matrix{Float64}(undef, length(residuals) - ar_len, ar_len)
for i in 1:ar_len
residuals_mat[:, i] = residuals[ar_len-i+1:n-i]
end
# Solve autoregressive problem
ar_coef = Variable(ar_len)
problem =
minimize(sumsquares(residuals_mat * ar_coef - residuals[ar_len+1:end]))
solve!(
problem,
MOI.OptimizerWithAttributes(ECOS.Optimizer, "maxit" => 200, "verbose" => 0),
)
```

```
Problem statistics
problem is DCP : true
number of variables : 1 (5 scalar elements)
number of constraints : 0 (0 scalar elements)
number of coefficients : 21_859
number of atoms : 3
Solution summary
termination status : NUMERICAL_ERROR
primal status : OTHER_RESULT_STATUS
dual status : OTHER_RESULT_STATUS
Expression graph
minimize
└─ qol (convex; positive)
├─ + (affine; real)
│ ├─ * (affine; real)
│ │ ├─ …
│ │ └─ …
│ └─ 3643×1 Matrix{Float64}
└─ [1;;]
```

Plot autoregressive fit of daily fluctuations for a few days:

```
ar_range = 1:145
day_range = ar_range .+ ar_len
plot(
day_range,
residuals[day_range],
label = "fluctuations from smooth fit",
ylabel = "Temperature difference (°C)",
)
plot!(
day_range,
residuals_mat[ar_range, :] * evaluate(ar_coef),
label = "autoregressive estimate",
xlabel = "Time (days)",
)
```

Now, we can add our autoregressive model for the residual temperatures to our smooth model to get an better fitting model for the daily temperatures in the city of Melbourne:

```
total_estimate = evaluate(yearly)
total_estimate[ar_len+1:end] += residuals_mat * evaluate(ar_coef)
```

```
3643-element Vector{Float64}:
15.663902485763881
15.344761935602632
15.354001303670024
16.134089753263478
18.438023944823033
17.019598399706254
15.289967968678694
14.188537906222022
16.39305927572535
18.46930145484137
⋮
13.84145223953454
14.098373964752284
11.868769604896553
13.873343052384929
14.409261498013906
13.955306048471076
13.796302824692958
13.904183045303483
15.167518253002223
```

We can plot the final fit of data across the whole time range:

```
plot(1:n, temps, label = "data", ylabel = "Temperature (°C)")
plot!(
1:n,
total_estimate,
label = "estimate",
xticks = 0:365:n,
xlabel = "Time (days)",
)
```

The RMS error of this final model is $\sim 2.3$:

```
root_mean_square_error =
sqrt(sum(x -> x^2, total_estimate - temps) / length(temps))
```

`2.3587029392911454`

*This page was generated using Literate.jl.*