Stationary Time Series

First of all, we need to define what is Stationary Time Series and figure out why we study Stationary Time Series.  Even before that, what is time series? According to the textbook, a time series is a set of obserations (Subscript[x, t]), each one being recorded at a specific time t. There are two types of Time Series. One is a discrete time series and the other is continuous time series. If the observations are over a time interval that is fixed, it’s discrete and if the observations are recorded continously over some time interval, it’s continuous.
Now, definition 1.3.1 from the textbook says that a time series model for the observed data {(Subscript[x, 1])} is a specification of the joint distributions (or possibly only the means and covariances) of a sequence of random variables {(Subscript[x, t])}of which {(Subscript[x, t])} is postulated to be a realization.

To state if the time series is stationary or not, we need to define the mean function and the covariance function first.
Definition 1.4.1
Let {(Subscript[X, t])} be a time series with E((Subscript[X, t]^2))<∞. The mean function of {(Subscript[X, t])} is
(Subscript[μ, X])(t) = E((Subscript[X, t]))

The covariance function of {(Subscript[X, t])} is 
(Subscript[γ, X])(r,s)= Cov((Subscript[X, r]),(Subscript[X, s]))=E[((Subscript[X, r])-(Subscript[μ, X])(r))((Subscript[X, s])-(Subscript[μ, X])(s))]
for all integers r and s.

To figure out if the time series is stationary or not, we need to use definition 1.4.2 saying 
{(Subscript[X, t])} is (weakly) stationary if 
(i) (Subscript[μ, X])(t) = E((Subscript[X, t])).
and 
(ii) (Subscript[γ, X])(t+h,t) is independent of t for each h.

Noise

(* sample time series of a roll of a dice with 6 faces. Rolled 80 times *)
data=RandomVariate[DiscreteUniformDistribution[{1,6}],80];
ListPlot[data,Filling->Axis]

(* Dice rolled many more times, the horizontal bands are exactly 1, 2, 3, 4, 5, 6 *)
data=RandomVariate[DiscreteUniformDistribution[{1,6}],8000];
ListPlot[data,Filling->Axis]

Now, Calculate the mean and variance of the distribution.

{Mean[data],Variance[data]}//N

h is the sliding window of the model that determines whether the stationary definition of the time series holds or not.
the code Length[data]-h and

t = 30;
(*timeshift by h *)
h = 1000;

datatimet = data[[t;;Length[data]-h]];

datatimetplush = data[[t+h;;Length[data]]];
{Mean[datatimet],Covariance[datatimet,datatimetplush]}//N

(* change this time see if the covariance changes *)
t = 190;

datatimet = data[[t;;Length[data]-h]];

datatimetplush = data[[t+h;;Length[data]]];
{Mean[datatimet],Covariance[datatimet,datatimetplush]}//N

t = 590;

datatimet = data[[t;;Length[data]-h]];

datatimetplush = data[[t+h;;Length[data]]];
{Mean[datatimet],Covariance[datatimet,datatimetplush]}//N

Let’s use Poisson instead of constant distribution for the dice

data=RandomVariate[PoissonDistribution[10],100];
ListPlot[data,Filling->Axis]

(* Increase the number of samples in time series *)
data=RandomVariate[PoissonDistribution[10],10^4];
ListPlot[data,Filling->Axis]

We do the same procedure as above. Calculate the mean and the variance of the data and vary t as we give any number of h to see if the sample of  time series is stationary or not. As we change the time, we can observe whether the covariance changes to check if the series is stationary.

{Mean[data],Variance[data]}//N

t = 30;
(*timeshift by h *)
h = 1000;

datatimet = data[[t;;Length[data]-h]];

datatimetplush = data[[t+h;;Length[data]]];
{Mean[datatimet],Covariance[datatimet,datatimetplush]}//N

(* change this time see if the covariance changes *)
t = 190;

datatimet = data[[t;;Length[data]-h]];

datatimetplush = data[[t+h;;Length[data]]];
{Mean[datatimet],Covariance[datatimet,datatimetplush]}//N

t = 590;

datatimet = data[[t;;Length[data]-h]];

datatimetplush = data[[t+h;;Length[data]]];
{Mean[datatimet],Covariance[datatimet,datatimetplush]}//N

Now let’s use another time series which is obtained frmo Normal distrubution, instead of rolling dice. Mean = 1 and Var = 3 
We will have one plot that has small sample and another with large sample to see the pattern.

data=RandomVariate[NormalDistribution[1,3],100];
ListPlot[data,Filling->Axis]

(* Now choose much large sample for time series *)
data=RandomVariate[NormalDistribution[1,3],10^4];
ListPlot[data,Filling->Axis]

Once again, we repeat the procedure as above, calculate mean and variance and the covariance.

{Mean[data],Variance[data]}//N

t= 900;
(*timeshift by h *)
h = 1000;

datatimet = data[[t;;Length[data]-h]];

datatimetplush = data[[t+h;;Length[data]]];
{Mean[datatimet],Covariance[datatimet,datatimetplush]}//N

(* change this time see if the covariance changes *)
t = 190;

datatimet = data[[t;;Length[data]-h]];

datatimetplush = data[[t+h;;Length[data]]];
{Mean[datatimet],Covariance[datatimet,datatimetplush]}//N

t = 590;

datatimet = data[[t;;Length[data]-h]];

datatimetplush = data[[t+h;;Length[data]]];
{Mean[datatimet],Covariance[datatimet,datatimetplush]}//N

All of the samples above had reasonably small changes to their mean and covariance as the time changed which we can conclude that all of them were (weakly) stationary models.

Random Walk

To study the random walk, it’s a good idea to cover the textbook definition first. 
Example 1.3.3 Random Walk
The random walk {(Subscript[S, t]),t=0,1,2,...}(starting at zero) is obtained by cumulatively summing (or “integrating”) iid random variables. Thus a random walk with zero mean is obtained by defining (Subscript[S, 0])=0 and (Subscript[S, t])= (Subscript[X, 1])+(Subscript[X, 2])+...+(Subscript[X, t]), for t=1,2,... where {(Subscript[X, t])} is iid noise. If  {(Subscript[X, t])}  is the binary process of example 1.3.2, then {(Subscript[S, t]),t=0,1,2,...} is called a simple symmetric random walk. This walk can be viewed as the location of a pedestrian who starts at position zero at time zero and at each integer time tosses a fair coin, stepping one unit to the right each time a head appears, and one unit to the left for each tail. A realization of length 200 of a simple symmetric random walk is shown in Figure 1.7. Notice that the outcomes of the coin tosses can be recovered from {(Subscript[S, t]),t=0,1,2,...} by differencing, THus the result of the tth toss can be found from (Subscript[S, t])-(Subscript[S, t-1])=(Subscript[X, t]).

this is the example that covers the random walk.
Example 1.4.3 The Random Walk
If {(Subscript[S, t])}is the random walk defined in Example 1.3.3 with {(Subscript[X, t])} as in Example 1.4.1 then (Subscript[ES, t])= 0, E((Subscript[S, t]^2))=((tσ)^2)<∞ for all t, and for h≥0,
(Subscript[γ, t])(t+h,t)=Cov((Subscript[S, t+h]),(Subscript[S, t]))
                       = Cov((Subscript[S, t])+(Subscript[X, t+1])+...+(Subscript[X, t+h,])(Subscript[S, t]))
                       = Cov((Subscript[S, t])+(Subscript[S, t]))
                       = ((tσ)^2).
Since (Subscript[γ, s])(t+h,t) depends on t, series {(Subscript[S, t])}is not stationary.

We now start with our own example. The model starts walking with flip a coin and H means go right +1 and T means go left -1. Take a step each time after the toss. 
Lets start with a 50 tosses.

(* start walking but at each step flip a coin, H go right +1, T go left -1, *)
data = RandomChoice[{-1,1},50];
dataaccum = Accumulate[data];
ListPlot[data]
ListLinePlot[dataaccum]

This graph gives you the distribustion of the accumlative steps you took from the origin:

Now that we have some distribution in 50 trials, lets try the large sample, 5000.

data = RandomChoice[{-1,1},5000];
dataaccum = Accumulate[data];
ListPlot[data]
ListLinePlot[dataaccum]

We can see some trends on the graph above. Now we calculate the mean and covariance with fixed h in different times.

{Mean[data],Variance[data]}//N

data = dataaccum;
t = 30;
(*timeshift by h *)
h = 150;

datatimet = data[[t;;Length[data]-h]];

datatimetplush = data[[t+h;;Length[data]]];
{Mean[datatimet],Covariance[datatimet,datatimetplush]}//N

(* change this time see if the covariance changes, as you can see it does!!! So random walk is NOT stationary *)
t = 90;

datatimet = data[[t;;Length[data]-h]];

datatimetplush = data[[t+h;;Length[data]]];
{Mean[datatimet],Covariance[datatimet,datatimetplush]}//N

t = 190;

datatimet = data[[t;;Length[data]-h]];

datatimetplush = data[[t+h;;Length[data]]];
{Mean[datatimet],Covariance[datatimet,datatimetplush]}//N

The result shows that the mean and variance vary as the time changes. Therefore, Random Walk is not a stationary model.

Auto Correlation

Autocorrelation function is the function we can use to help finding the covariance of the time series. Here is the definition.
Definition 1.4.3 The autocorrelation function (ACF) of {(Subscript[X, t])} at lag h is
(Subscript[ρ, X])(h)=(((Subscript[γ, X])(h))/((Subscript[γ, X])(0)))= Cor((Subscript[X, t+h]),(Subscript[X, t])).
Also, the definition 1.4.4 gives the definition of the sample correlation function
(OverHat[ρ])(h)=(((OverHat[γ])(h))/((OverHat[γ])(0))), -n<h<n.

The autocorrealtion function has to be constant in order to be independent of t. If it is independent of t, the series is stationary.

So we assign a function naming “autoCorrelation” as below.

(* Auto Correlation (Subscript[ρ, X])(h). However the time t > 0 is added for the sake of emphasis and in case the time series is not stationary.
This function is actually called SAMPLE Auto Correlation *)
autoCorrelation[data0_, h0_, t0_]:= Module[{data = data0, h =h0, t =t0,datatimet,datatimetplush,cov1h, cov10},

datatimet = data[[t;;Length[data]-h]];

datatimetplush = data[[t+h;;Length[data]]];
cov1h=Covariance[datatimet,datatimetplush];

h=0;
datatimet = data[[t;;Length[data]-h]];
datatimetplush = data[[t+h;;Length[data]]];
cov10=Covariance[datatimet,datatimetplush];

cov1h/cov10//N
]

We use the dice rolling example used above to see the correlation works on identifying the stationariness.

(* Dice rolled many more times, the horizontal bands are exactly 1, 2, 3, 4, 5, 6 *)
data=RandomVariate[DiscreteUniformDistribution[{1,6}],8000];
{Mean[data],Variance[data]}//N

autoCorrelation[data, 900, 600]

h = 1000;
ListPlot[Table[autoCorrelation[data, h,t], {t, 1, Ceiling[Length[data]/2]}], PlotRange->{{0,4000}, {0,1}}]

Both the value of the autocorrelation and the plot shows that it’s constant. Therefore, the series is stationary.

We use the poisson distribution once more.

data=RandomVariate[PoissonDistribution[10],10^4];
ListPlot[data,Filling->Axis]

h = 1000;
ListPlot[Table[autoCorrelation[data, h,t], {t, 1, Ceiling[Length[data]/2]}], PlotRange->{{0,4000}, {0,1}}]

the autocorrelation is constant again, thus stationary.


Create another random walk

data = RandomChoice[{-1,1},3000];
dataaccum = Accumulate[data];
ListPlot[data]
ListLinePlot[dataaccum]

As you can see its Auto Correlation function is not constant in time t

h = 1000;
ListPlot[Table[autoCorrelation[dataaccum, h,t], {t, 1, Ceiling[Length[dataaccum]/2]}], PlotRange->{{0,Length[dataaccum]}, {-1,1}}]

Classical Decomposition Model

Section 1.5 emphasize on the estimation and elimination of trend and seasonal components.
First to study is a classical decomposition model.
Classical Decomposition Model
(Subscript[X, t])= (Subscript[m, t])+(Subscript[s, t])+(Subscript[Y, t]), t=1,2,...,n,
where (Subscript[EY, t])=0, (Subscript[s, t+d])= (Subscript[s, t]), and (∑_(j=1)^d)(Subscript[s, j])=0.
The below plots are some of the example models with 1. slowly changing trend, seasonal component (periodic), and random noise. and the last one is the composition of all of those functions.

t = .;
m=.;
s=.;
X=.;
Y=.;
(* Trend, slowly changing function *)
m = 0.1*t^2;
(* Seasonal Component, preodic cyclic *)
s = 7*Cos[t];
(* Random Noise *)
Y = HoldForm[4*RandomReal[{-1,1}]];

X[w_]:= (m+s+ReleaseHold[Y])/.{t->w};

Plot[m, {t, 0, 50}]
Plot[s, {t, 0, 50}]
Plot[ReleaseHold[Y], {t, 0, 50}]
Plot[X[t], {t, 0, 50}]

Smoothing

Sometimes, we need to smooth the time series to get a better sense of the trend because some people might get confused with the noise and the trend. To reduce the chance of error, lets study some of the examples of the smoothing. Firt of all, we assign a Xseries which is the sequence with difference of 0.3.

(* Grab a signal and set the time sequence different at 0.3.
Look at the data to get an idea. *)

Xseries = Table[X[t], {t, 0,50, 0.3}]
ListPlot[Xseries]

Moving Average

The first method to cover is the moving average. This is the definition from the book.
Smoothing with a finite moving average filter.
Let q be a nonnegative integer and consider the two-sided moving average
(Subscript[W, t])=(((2q+1))^(-1))(∑_(j=-q)^q)(Subscript[X, t-j])
of the process {(Subscript[X, t])} defined by (Subscript[X, t])= (Subscript[m, t])+(Subscript[Y, t]), t=1,...,n, where (Subscript[EY, t])=0.  Then for q+1≤ t≤n-q,
(Subscript[W, t])=(((2q+1))^(-1))(∑_(j=-q)^q)(Subscript[m, t-j])+(((2q+1))^(-1))(∑_(j=-q)^q)(Subscript[Y, t-j])≈(Subscript[m, t]),
assuming that (Subscript[m, t]) is approximately linear over the interval [t-q,t+q] and that the average of the error terms over this interval is close to zero.
The moving average thus provides us with the estimates
(Subscript[(OverHat[m]), t])=(((2q+1))^(-1))(∑_(j=-q)^q)(Subscript[X, t-j]),q+1≤t≤n-q.

So we perform the moving average below to see if there is any change to the plot.

MovingAverage[{a,b,c,d,e,f,g}, 3]

movingaverage=MovingAverage[Xseries, 4];
ListPlot[movingaverage]

The plot flows better than the original graph.

Exponential Moving Average (Smoothing) 



ExponentialMovingAverage[x,α] generates a list of results y in which (Subscript[y, i+1])=(Subscript[y, i])+α((Subscript[x, i+1])-(Subscript[y, i])).

Let study the definition of the exponential smoothing. 
For any fixed α ∈[0,1], the one sided moving average (Subscript[(OverHat[m]), t]), t=1,...,n, defined by the recursions


(Subscript[(OverHat[m]), t]) = (Subscript[αX, t])+(1-α)(Subscript[(OverHat[m]), t-1]), t=2,...,n, and (Subscript[(OverHat[m]), 1])= (Subscript[X, 1])

ExponentialMovingAverage[{a,b,c},α]

movingaverage=ExponentialMovingAverage[Xseries, .2];
ListPlot[movingaverage]

Differencing


If list has length s, Differences[list,n] has length s-n.

We can also eliminate the trend with the method of differencing.
According to the textbook,
Instead of attempting to remove the noise by smoothing we now attempt to eliminate the trend by differencing. 
So this is a different method from smoothing techniques as above. We define the lag-1 difference operator  ▽ by 
(Subscript[▽X, t])= (Subscript[X, t])-(Subscript[X, t-1])=(1-B)(Subscript[X, t]), where B is the backward shift operator,
(Subscript[BX, t])=(Subscript[X, t-1]).

Differences[{a,b,c,d,e}]

Differences[{a,b,c,d,e},2]

power=5
diff = Differences[Xseries, power];
ListPlot[Drop[Xseries,-power]]

By applying the the differencing method, the graph is detrended even though it’s not as smooth as the smoothing techniques.