Linear Stochastic Models

Special Types of Random Processes: AR,

MA, and ARMA

Digital Signal Processing

Department of Electrical and Electronic Engineering, Imperial College

[email protected]

c Danilo P. Mandic
Digital Signal Processing 1
Motivation:- Wold Decomposition Theorem

The most fundamental justification for time series analysis is due to Wold’s
decomposition theorem, where it is explicitly proved that any (stationary)
time series can be decomposed into two different parts.
Therefore, a general random process can be written a sum of two processes

x[n] = xp[n] + xr [n]

⇒ xr [n] – regular random process

⇒ xp[n] – predictable process, with xr [n] ⊥ xp[n],

E{xr [m]xp[n]} = 0

that is we can separately treat the predictable process (i.e. a

deterministic signal) and a random signal.

c Danilo P. Mandic
Digital Signal Processing 2
What do we actually mean?

a) Periodic oscillations b) Small nonlinearity c) Route to chaos

d) Route to chaos e) small noise f) HMM and others

c Danilo P. Mandic
Digital Signal Processing 3
Example from brain science

Electrode positions Raw EEG Useful signal

c Danilo P. Mandic
Digital Signal Processing 4
Linear Stochastic Processes
It therefore follows that the general form for the power spectrum of a WSS
process is

N
X
Px(eω ) = Pxr (eω ) + αk u0(ω − ωk )
k=1

We look at processes generated by filtering white noise with a linear

shift–invariant filter that has a rational system function. These include the

• Autoregressive (AR) → all pole system

• Moving Average (MA) → all zero system

• Autoregressive Moving Average (ARMA) → poles and zeros

Notice the difference between shift–invariance and time–invariance

c Danilo P. Mandic
Digital Signal Processing 5
 ACF and Spectrum of ARMA models
Much of interest are the autocorrelation function and power spectrum of
these processes. (Recall that ACF ≡ P SD in terms of the available
information)
Suppose that we filter white noise w[n] with a causal linear shift–invariant
filter having a rational system function with p poles and q zeros
Pq
Bq (z) k=0 bq (k)z −k
H(z) = = Pp
Ap(z) 1 + k=1 ap(k)z −k

Assuming that the filter is stable, the output process x[n] will be
2
wide–sense stationary and with Pw = σw , the power spectrum of x[n] will
be
−1
2 Bq (z)Bq (z )
Px(z) = σw
Ap(z)Ap(z −1)
∗
Recall that “(·) ” in analogue frequency corresponds to “z −1” in “digital freq.”

c Danilo P. Mandic
Digital Signal Processing 6
Frequency Domain
In terms of “digital” frequency θ (unit circle – e−θ = e−ωT )

• Bq (z)Bq (z −1) # “quadratic form” and real valued

• Ap(z)Ap(z −1) # “quadratic form” and real valued

Bq (eθ )2
 
Pz (eθ ) = σw
2
2
|Ap(eθ )|

We are therefore using H(z) to shape the spectrum of white noise.

A process having a power spectrum of this form is known as an
autoregressive moving average process of order (p, q) and is referred to
as an
ARMA(p,q) process

c Danilo P. Mandic
Digital Signal Processing 7
Example
Plot the power spectrum of an ARMA(2,2) process for which
• the zeros of H(z) are z = 0.95e±π/2

• poles are at z = 0.9e±2π/5

Solution: The system function is (poles and zeros – resonance & sink)
1 + 0.9025z −2
H(z) =
1 − 0.5562z −1 + 0.81z −2
7
Power Spectrum

−1
0 0.5 1 1.5 2 2.5 3 3.5
Frequency

c Danilo P. Mandic
Digital Signal Processing 8
Difference Equation Representation
Random processes x[n] and w[n] are related by the linear constant
coefficient equation

p
X q
X
x[n] − ap(l)x[n − l] = bq (l)w[n − l]
l=1 l=0

Notice that the autocorrelation function of x[n] and crosscorrelation

between x[n] and w[n] follow the same difference equation, i.e. if we
multiply both sides of the above equation by x[n − k] and take the
expected value, we have

p
X q
X
rxx(k) − ap(l)rxx(k − l) = bq (l)rxw (k − l)
l=1 l=0

Since x is WSS, it follows that x[n] and w[n] are jointly WSS.

c Danilo P. Mandic
Digital Signal Processing 9
General Linear Processes: Stationarity and Invertibility
Consider a linear stochastic process # output from a linear filter, driven by
WGN w[n]
∞
X
x[n] = w[n] + b1w[n − 1] + b2w[n − 2] + · · · = w[n] + bj w[n − j]
j=1

that is, a weighted sum of past inputs w[n].

For this process to be a valid P
stationary process, the coefficients must be
absolutely summable, that is ∞ j=0 |bj | < ∞.

The model implies that under suitable condition, x[n] is also a weighted
sum of past values of x, plus an added shock w[n], that is

x[n] = a1x[n − 1] + a2x[n − 2] + · · · + w[n]

P∞
• Linear Process is stationary if j=0 |bj | < ∞
P∞
• Linear Process is invertible if j=0 |aj | < ∞

c Danilo P. Mandic
Digital Signal Processing 10
Are these ARMA(p,q) processes?


0, n < 0
• Unit response u[n] =
1, n ≥ 0
– If w[n] = δ[n] then

u[n] = u[n − 1] + w[n], n≥0


0, n < 0
• Ramp function r[n] =
n, n ≥ 0
– If w[n] = u[n] then

r[n] = r[n − 1] + w[n], n≥0

c Danilo P. Mandic
Digital Signal Processing 11
Autoregressive Processes

A general AR(p) process (autoregressive of order p) is given by

p
X
x[n] = a1x[n − 1] + · · · + apx[n − p] + w[n] = aix[n − i] + w[n]
i=1
Observe the auto–regression above

Duality between AR and MA processes:

For instance the first order autoregressive process
∞
X
x[n] = a1x[n − 1] + w[n] ⇔ bj w[n − j]
j=0

Due to its “all–pole“ nature follows the duality between IIR and FIR filters.

c Danilo P. Mandic
Digital Signal Processing 12
ACF and Spectrum of AR Processes
To obtain the autocorrelation function of an AR process, multiply the
above equation by x[n − k] to obtain

x[n − k]x[n] = a1x[n − k]x[n − 1] + a2x[n − k]x[n − 2] + · · ·

+apx[n − k]x[n − p] + x[n − k]w[n]

Notice that E{x[n − k]w[n]} vanishes when k > 0. Therefore we have

rxx(k) = a1rxx(k − 1) + a2rxx(k − 2) + · · · + aprxx(k − p) k > 0

On dividing throughout by rxx(0) we obtain

ρ(k) = a1ρ(k − 1) + a2ρ(k − 2) + · · · + apρ(k − p) k > 0

Parameters ρ(k) are correlation coefficients

c Danilo P. Mandic
Digital Signal Processing 13
Variance and Spectrum of AR Processes
Variance:
2
When k = 0 the contribution from the term E{x[n − k]w[n]} is σw , and

2
rxx(0) = a1rxx(−1) + a2rxx(−2) + · · · + aprxx(−p) + σw

Divide by rxx(0) = σx2 to obtain

2
σ w
σx2 =
1 − ρ1a1 − ρ2a2 − · · · − ρpap

Spectrum:

2
2σw
Pxx(f ) = 2 0 ≤ f ≤ 1/2
|1 − a1 e−2πf − · · · − ap e−2πpf |

Recall Spectrum of linear systems from the Course Introduction

c Danilo P. Mandic
Digital Signal Processing 14
 Yule–Walker Equations
For k = 1, 2, . . . , p from the general autocorrelation function, we obtain
a set of equations:-

rxx(1) = a1rxx(0) + a2rxx(1) + · · · + aprxx(p − 1)

rxx(2) = a1rxx(1) + a2rxx(0) + · · · + aprxx(p − 2)
.. = ..

rxx(p) = a1rxx(p − 1) + a2rxx(p − 2) + · · · + aprxx(0)

These equations are called the Yule–Walker or normal equations.

Their solution gives us the set of autoregressive parameters
a = [a1, . . . , ap]T . This can be expressed in a vector–matrix form as

a = R−1
xx rxx

Due to Toeplitz structure of Rxx, its positive definitness enables matrix inversion

c Danilo P. Mandic
Digital Signal Processing 15
ACF Coefficients

For the autocorrelation coefficients

ρk = rxx(k)/rxx(0)

we have

ρ1 = a1 + a2ρ1 + · · · + apρp−1
ρ2 = a1ρ1 + a2 + · · · + apρp−2
.. = ..

ρp = a1ρp−1 + a2ρp−2 + · · · + ap

When does the sequence {ρ0, ρ1, ρ2, . . .} vanish?

Homework:- Try command xcorr in Matlab

c Danilo P. Mandic
Digital Signal Processing 16
Example:- Yule–Walker modelling in Matlab
In Matlab – Power spectral density using Y–W method pyulear

Pxx = pyulear(x,p)
[Pxx,w] = pyulear(x,p,nfft)
[Pxx,f] = pyulear(x,p,nfft,fs)
[Pxx,f] = pyulear(x,p,nfft,fs,’range’)
[Pxx,w] = pyulear(x,p,nfft,’range’)

Description:-

Pxx = pyulear(x,p)

implements the Yule-Walker algorithm, and returns Pxx, an estimate of the

power spectral density (PSD) of the vector x.

To remember for later → This estimate is also an estimate of the

maximum entropy.
Se also aryule, lpc, pburg, pcov, peig, periodogram

c Danilo P. Mandic
Digital Signal Processing 17
Example:- AR(p) signal generation
• Generate the input signal x by filtering white noise through the AR
filter

• Estimate the PSD of x based on a fourth-order AR model

Solution:-
randn(’state’,1);
x = filter(1,a,randn(256,1)); % AR system output
pyulear(x,4) % Fourth-order estimate

c Danilo P. Mandic
Digital Signal Processing 18
Alternatively:- Yule–Walker modelling
AR(4) system given by
y[n] = 2.2137y[n−1]−2.9403y[n−2]+2.1697y[n−3]−0.9606y[n−4]+w[n]
a = [1 -2.2137 2.9403 -2.1697 0.9606]; % AR filter coefficients
freqz(1,a) % AR filter frequency response
title(’AR System Frequency Response’)

c Danilo P. Mandic
Digital Signal Processing 19
From Data to AR(p) Model

So far, we assumed the model (AR, MA, or ARMA) and analysed the ACF
and PSD based on known model coefficients.
In practice:- DATA # MODEL

This procedure is as follows:-

* record data x(k)

* find the autocorrelation of the data ACF(x)
* divide by r_xx(0) to obtain correlation coefficients \rho(k)
* write down Yule-Walker equations
* solve for the vector of AR paramters

The problem is that we do not know the model order p beforehand;

we will deal with this problem later in Lecture 2.

c Danilo P. Mandic
Digital Signal Processing 20
Example:- Finding parameters of
x[n] = 1.2x[n − 1] − 0.8x[n − 2] + w[n]

AR(2) signal x=filter([1],[1, −1.2, 0.8],w) ACF for AR(2) signal x=filter([1],[1, −1.2, 0.8],w) ACF for AR(2) signal x=filter([1],[1, −1.2, 0.8],w)
6 2000

1500
1500
4

ACF of AR(2) signal

ACF of AR(2) signal

AR(2) signal values

1000
1000
2

500 500

0
0 0

−2
−500
−500

−4
−1000
−1000

−6 −1500
0 50 100 150 200 250 300 350 400 −400 −300 −200 −100 0 100 200 300 400 −20 −10 0 10 20
Sample number Correlation lag Correlation lag

Apply:- for i=1:6; [a,e]=aryule(x,i); display(a);end

a(1) = [0.6689] a(2) = [1.2046, −0.8008]

a(3) = [1.1759, −0.7576, −0.0358]
a(4) = [1.1762, −0.7513, −0.0456, 0.0083]
a(5) = [1.1763, −0.7520, −0.0562, 0.0248, −0.0140]
a(6) = [1.1762, −0.7518, −0.0565, 0.0198, −0.0062, −0.0067]

c Danilo P. Mandic
Digital Signal Processing 21
 Special case:- AR(1) Process (Markov)
Given below (Recall p(x[n], x[n − 1], . . . , x[0]) = p(x[n] | x[n − 1]))

x[n] = a1x[n − 1] + w[n] = w[n] + a1x[n − 1] + a21w[n − 2] + · · ·

i) for the process to be stationary −1 < a1 < 1.

ii) Autocorrelation Function:- from Yule-Walker equations

rxx(k) = a1rxx(k − 1), k>0

or for the correlation coefficients, with ρ0 = 1

ρk = ak1 , k>0

Notice the difference in the behaviour of the ACF for a1 positive and negative

c Danilo P. Mandic
Digital Signal Processing 22
Variance and Spectrum of AR(1) process
Can be calculated directly from a general expression of the variance and
spectrum of AR(p) processes.

• Variance:- Also from a general expression for the variance of linear

processes from Lecture 1

2 2
σw σw
σx2 = =
1 − ρ1a1 1 − a21

• Spectrum:- Notice how the flat PSD of WGN is shaped according to

the position of the pole of AR(1) model (LP or HP)

2 2
2σw 2σw
Pxx(f ) = 2 = 1 + a2 − 2a cos(2πf )
|1 − a1e−2πf | 1
1

c Danilo P. Mandic
Digital Signal Processing 23
Example: ACF and Spectrum of AR(1) for a = ±0.8
x[n] = −0.8*x[n−1] + w[n] x[n] = 0.8*x[n−1] + w[n]
4 5

Signal values

Signal values
0 0

−2

−4 −5
0 20 40 60 80 100 0 20 40 60 80 100
Sample Number Sample Number
ACF ACF
1 1

0.5
Correlation

Correlation
0 0.5

−0.5

−1 0
0 5 10 15 20 0 5 10 15 20
Correlation lag Correlation lag
Power/frequency (dB/rad/sample)

Power/frequency (dB/rad/sample)
Burg Power Spectral Density Estimate Burg Power Spectral Density Estimate
10 15

5 10

5
0
0
−5
−5
−10 −10

−15 −15
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1
Normalized Frequency (×π rad/sample) Normalized Frequency (×π rad/sample)

a < 0 → High Pass a > 0 → Low Pass

c Danilo P. Mandic
Digital Signal Processing 24
Special Case:- Second Order Autoregressive Processes
AR(2)

The input–output functional relationship is given by

x[n] = a1x[n − 1] + a2x[n − 2] + w[n]

For stationarity- (to be proven later)

a1 + a2 < 1
a2 − a1 < 1
−1 < a2 < 1

This will be shown within the so–called “stability triangle”

c Danilo P. Mandic
Digital Signal Processing 25
Work by Yule – Modelling of sunspot numbers

Recorded for more than 300 years.

In 1927, Yule modelled them and invented AR(2) model
Sunspot series
150

100

Signal values
50

−50
0 50 100 150 200 250 300
Sample Number

ACF for sunspot series

1

0.5
Correlation

−0.5
0 5 10 15 20 25 30 35 40 45 50
Correlation lag

Sunspot numbers and its autocorrelation function

c Danilo P. Mandic
Digital Signal Processing 26
Autocorrelation function of AR(2) processes
The ACF
ρk = a1ρk−1 + a2ρk−2 k>0

• Real roots: ⇒ (a21 + 4a2 > 0) ACF = mixture of damped exponentials

• Complex roots: ⇒ (a21 + 4a2 < 0) ⇒ ACF exhibits a pseudo–periodic
behaviour
D k sin(2πf0 k + F )
ρk =
sin F
D - damping factor, of a sine wave with frequency f0 and phase F.
√
D = −a2
a1
cos(2πf0) = √
2 −a2
1 + D2
tan(F ) = tan(2πf0 )
1 − D2

c Danilo P. Mandic
Digital Signal Processing 27
Stability Triangle
a2

1
ACF ACF

II I

m m
Real Roots

−2 2 a1
ACF
ACF

m
m
III Complex Roots IV

−1

i) Real roots Region 1: Monotonically decaying ACF

ii) Real roots Region 2: Decaying oscillating ACF
iii) Complex roots Region 3: Oscilating pseudoperiodic ACF
iv) Complex roots Region 4: Pseudoperiodic ACF

c Danilo P. Mandic
Digital Signal Processing 28
Yule–Walker Equations
Substituting p = 2 into Y-W equations we have
ρ1 = a1 + a2ρ1
ρ2 = a1ρ1 + a2

which when solved for a1 and a2 gives

ρ1(1 − ρ2) ρ2 − ρ21
a1 = a2 =
1 − ρ21 1 − ρ21
or substituting in the equation for ρ

a1
ρ1 =
1 − a2
a21
ρ2 = a2 +
1 − a2

c Danilo P. Mandic
Digital Signal Processing 29
Variance and Spectrum

More specifically, for the AR(2) process, we have:-

Variance
2 2
 
2 σw 1 − a2 σw
σx = =
1 − ρ1a1 − ρ2a2 1 + a2 (1 − a2)2 − a21

Spectrum

2
2σw
Pxx(f ) = 2
|1 − a1e−2πf − a2e−4πf |
2
2σw
= , 0 ≤ f ≤ 1/2
1 + a21 + a22 − 2a1(1 − a2 cos(2πf ) − 2a2 cos(4πf ))

c Danilo P. Mandic
Digital Signal Processing 30
Example AR(2): x[n] = 0.75x[n − 1] − 0.5x[n − 2] + w[n]

x[n] = 0.75*x[n−2] − 0.5*x[n−1] + w[n]

20

10

Signal values
0

−10

−20
0 50 100 150 200 250 300 350 400 450 500
Sample Number
ACF
1
Correlation

0.5

−0.5
0 5 10 15 20 25 30 35 40 45 50
Correlation lag
Power/frequency (dB/rad/sample)

Burg Power Spectral Density Estimate

15

10

−5

−10
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

√ Normalized Frequency (×π rad/sample)

cos−1 (0.5303) 1
The damping factor D = 0.5 = 0.71, frequency f0 = 2π = 6.2
The fundamental period of the autocorrelation function is 6.2.

c Danilo P. Mandic
Digital Signal Processing 31
Partial Autocorrelation Function:- Motivation
Let us revisit example from page 21 of Lecture Slides.
AR(2) signal x=filter([1],[1, −1.2, 0.8],w) ACF for AR(2) signal x=filter([1],[1, −1.2, 0.8],w) ACF for AR(2) signal x=filter([1],[1, −1.2, 0.8],w)
6 2000

1500
1500
4

ACF of AR(2) signal

ACF of AR(2) signal

AR(2) signal values

1000
1000
2

500 500

0
0 0

−2
−500
−500

−4
−1000
−1000

−6 −1500
0 50 100 150 200 250 300 350 400 −400 −300 −200 −100 0 100 200 300 400 −20 −10 0 10 20
Sample number Correlation lag Correlation lag

We do not know p, let us re-write the coefficients as [a_1p,...,a_p

p = 1 # [0.6689] = a11 p = 2 # [1.2046, −0.8008] = [a21, a22]

p = 3 #[1.1759, −0.07576, −0.0358] = [a31, a32, a33]
p = 4 # [1.1762, −0.7513, −0.0456, 0.0083] = [a41, a42, a43, a44]
p = 5 # [1.1763, −0.7520, −0.0562, 0.0248, −0.0140] = [a51, . . . , a55]
p = 6 # [1.1762, −0.7518, −0.0565, 0.0198, −0.0062, −0.0067] =
[a61, . . . , a66]

c Danilo P. Mandic
Digital Signal Processing 32
Partial Autocorrelation Function

Notice: ACF of AR(p) infinite in duration, but can by be described in

terms of p nonzero functions ACFs.
Denote by akj the jth coefficient in an autoregressive representation of
order k, so that akk is the last coefficient. Then

ρj = akj ρj−1 + · · · + ak(k−1)ρj−k+1 + akk ρj−k j = 1, 2, . . . , k

leading to the Yule–Walker equation, which can be written as

    
1 ρ1 ρ2 ··· ρk−1 ak1 ρ1
 ρ1 1 ρ1 ··· ρk−2   ak2   ρ2 
= . 

 .. .. .. ... ..  
 .
 .   . 
ρk−1 ρk−2 ρk−3 ··· 1 akk ρk

c Danilo P. Mandic
Digital Signal Processing 33
Partial ACF Coefficients:

Solving these equations for k = 1, 2, . . . successively, we obtain

 

 1 ρ1 ρ1 


 ρ1 1 ρ2 

ρ2 − ρ21  ρ2 ρ1 ρ3 
a11 = ρ1, a22 = , a33 =  , etc
1 − ρ21  1 ρ1 ρ2 


 ρ1 1 ρ1 

 ρ2 ρ1 1 

• The quantity akk , regarded as a function of lag k, is called the partial

autocorrelation function.

• For an AR(p) process, the PAC akk will be nonzero for k ≤ p and zero
for k > p ⇒ tells us the order of an AR(p) process.

c Danilo P. Mandic
Digital Signal Processing 34
Importance of Partial ACF
For a zero mean process x[n], the best linear predictor in the mean
square error sense of x[n] based on x[n − 1], x[n − 2], . . . is
x̂[n] = ak−1,1x[n − 1] + ak−1,2x[n − 2] + · · · + ak−1,k−1 x[n − k + 1]

(apply the E{·} operator to the general AR(p) model expression, and
recall that E{w[n]} = 0)
(Hint:
E{x[n]} = x̂[n] = E {ak−1,1 x[n − 1] + · · · + ak−1,k−1 x[n − k + 1] + w[n]} =
ak−1,1 x[n − 1] + · · · + ak−1,k−1 x[n − k + 1]) )
whether the process is an AR or not

In MATLAB, check the function:

ARYULE
and functions
PYULEAR, ARMCOV, ARBURG, ARCOV, LPC, PRONY

c Danilo P. Mandic
Digital Signal Processing 35
Model order for Sunspot numbers

Sunspot series ACF for sunspot series

150 1

100

Signal values
0.5

Correlation
50

0
0

−50 −0.5
0 100 200 300 0 10 20 30 40 50
Sample Number Correlation lag

Partial ACF for sunspot series Burg Power Spectral Density Estimate

Power/frequency (dB/rad/sample)
1 35

30
0.5
25
Correlation

0 20

15
−0.5
10

−1 5
0 10 20 30 40 50 0 0.2 0.4 0.6 0.8 1
Correlation lag Normalized Frequency (×π rad/sample)

Sunspot numbers, their ACF and partial autocorrelation (PAC)

After lag k = 2, the PAC becomes very small

c Danilo P. Mandic
Digital Signal Processing 36
Model order for AR(2) generated process

AR(2) signal ACF for AR(2) signal

8 1

4 0.5

Signal values

Correlation
2

0 0

−2

−4 −0.5

−6

−8 −1
0 100 200 300 400 500 0 10 20 30 40 50
Sample Number Correlation lag

Partial ACF for AR(2) signal Burg Power Spectral Density Estimate

Power/frequency (dB/rad/sample)
0.8 15

0.6 10
0.4
5
Correlation

0.2
0
0
−5
−0.2
−10
−0.4

−0.6 −15

−0.8 −20
0 10 20 30 40 50 0 0.2 0.4 0.6 0.8 1
Correlation lag Normalized Frequency (×π rad/sample)

AR(2) signal, its ACF and partial autocorrelation (PAC)

After lag k = 2, the PAC becomes very small

c Danilo P. Mandic
Digital Signal Processing 37
Model order for AR(3) generated process

AR(3) signal ACF for AR(3) signal

10 1

5 0.8

Signal values

Correlation
0 0.6

−5 0.4

−10 0.2

−15 0
0 100 200 300 400 500 0 10 20 30 40 50
Sample Number Correlation lag

Partial ACF for AR(3) signal Burg Power Spectral Density Estimate

Power/frequency (dB/rad/sample)
0.4 30

0.2
20
0
Correlation

10
−0.2

−0.4
0

−0.6
−10
−0.8

−1 −20
0 10 20 30 40 50 0 0.2 0.4 0.6 0.8 1
Correlation lag Normalized Frequency (×π rad/sample)

AR(3) signal, its ACF and partial autocorrelation (PAC)

After lag k = 3, the PAC becomes very small

c Danilo P. Mandic
Digital Signal Processing 38
Model order for a financial time series
From:- http://finance.yahoo.com/q/ta?s=%5EIXIC&t=1d&l=on&z=m&q=b&p=v&a=&c=
Nasdaq ascending Nasdaq descending
Nasdaq composite June 2003 − February 2007 Nasdaq composite February 2007 − June 2003
2600 2600

2400 2400
Nasdaq value

Nasdaq value
2200 2200

2000 2000

1800 1800

1600 1600

1400 1400
0 500 1000 1500 2000 0 500 1000 1500 2000
Day number Day number
ACFx of
10 Nasdaq composite June 2003 − February 2007
7 ACFx of
10 Nasdaq composite June 2003 − February 2007
7

7 7

6 6

5 5

4 4
ACF value

ACF value
3 3

2 2

1 1

0 0

−1 −1

−2 −2

−3 −3
−2000 −1500 −1000 −500 0 500 1000 1500 2000 −2000 −1500 −1000 −500 0 500 1000 1500 2000
Correlation lag Correlation lag

c Danilo P. Mandic
Digital Signal Processing 39
Partial ACF for financial time series

a = 1.0000 -0.9994

a = 1.0000 -0.9982 -0.0011

a = 1.0000 -0.9982 0.0086 -0.0097

a = 1.0000 -0.9983 0.0086 -0.0128 0.0030

a = 1.0000 -0.9983 0.0086 -0.0128 0.0026 0.0005

a = 1.0000 -0.9983 0.0086 -0.0127 0.0026 0.0017 -0.0

c Danilo P. Mandic
Digital Signal Processing 40
Model Order Selection – Practical issues
In practice – the greater the model order the higher the accuracy
⇒ When do we stop?

To save on computational complexity, we introduce “penalty” for a high

model order. The criteria for model order selection are, for instance MDL
(minimum description length - Rissanen), AIC (Akaike Information
criterion), given by

p ∗ log(N )
M DL = log(E) +
N
AIC = log(E) + 2p/N

E = the loss function (typically cumulative squared error,

p = the number of estimated parameters
N = the number of estimated data.

c Danilo P. Mandic
Digital Signal Processing 41
Example:- Model order selection – MDL vs AIC
Let us have a look at the squared error and the MDL and AIC criteria for
an AR(2) model with
a1 = 0.5 a2 = −0.3
MDL for AR(2)
1
MDL
0.98 Cumulative Squared Error

0.96

0.94

0.92

0.9

0.88
1 2 3 4 5 6 7 8 9 10

AIC for AR(2)

1
AIC
0.98 Cumulative Squared Error

0.96

0.94

0.92

0.9

0.88
1 2 3 4 5 6 7 8 9 10
AR(2) Model Order

(Model error)2 versus the model order p

c Danilo P. Mandic
Digital Signal Processing 42
Moving Average Processes
A general MA(q) process is given by
x[n] = w[n] + b1w[n − 1] + · · · + bq w[n − q]

Autocorrelation function: The autocovariance function of MA(q)

ck = E[(w[n] + b1w[n − 1] + · · · + bq w[n − q])(w[n − k]
+b1w[n − k − 1] + · · · + bq w[n − k − q])]

Hence the variance of the process

c0 = (1 + b21 + · · · + b2q )σw
2

The ACF of an MA process has a cutoff after lag q.

Spectrum: All zeros ⇒ struggles to model PSD with peaks

2  −2πf −4πf −2πqf 
 
P (f ) = 2σw 1 + b1e + b2 e + · · · + bq e

c Danilo P. Mandic
Digital Signal Processing 43
Example:- MA(3) process

MA(3) signal ACF for MA(3) signal

3 1.2

1
2
0.8

Signal values

Correlation
1
0.6

0.4
0

0.2
−1
0

−2 −0.2
0 100 200 300 400 500 0 10 20 30 40 50
Sample Number Correlation lag

Partial ACF for MA(3) signal Burg Power Spectral Density Estimate

Power/frequency (dB/rad/sample)
0.5 −4

0.4
−6
0.3
−8
Correlation

0.2

0.1 −10

0
−12
−0.1
−14
−0.2

−0.3 −16
0 10 20 30 40 50 0 0.2 0.4 0.6 0.8 1
Correlation lag Normalized Frequency (×π rad/sample)

MA(3) model, its ACF and partial autocorrelation (PAC)

After lag k = 3, the ACF becomes very small

c Danilo P. Mandic
Digital Signal Processing 44
Analysis of Nonstationary Signals
Speech Signal
1
W1 W2 W3

Signal values
0.5

0
2000 3000 4000 5000 6000 7000 8000 9000 10000
Sample Number
Partial ACF for W1 Partial ACF for W2 Partial ACF for W3
0 1 1

−1 0.5 0.5
Correlation

Correlation

Correlation
1 0 0

0 −0.5 −0.5

−1 −1 −1
0 25 50 0 25 50 0 25 50
Correlation lag Correlation lag Correlation lag
MDL calculated for W1 MDL calculated for W2 MDL calculated for W3
1 1 1

0.8 Calculated Calculated Model 0.8 Calculated Model

Model Order = 13 Order > 50 Order = 24
MDL

MDL

MDL
0.6 0.5 0.6

0.4 0.4

0.2 0 0.2
0 25 50 0 25 50 0 25 50
Model Order Model Order Model Order

Different AR models for different segments of speech

To deal with nonstationarity we need short sliding windows

c Danilo P. Mandic
Digital Signal Processing 45
Duality Between AR and MA Processes

i) A stationary finite AR(p) process can be represented as an infinite order

MA process. A finite MA process can be represented as an infinite AR
process.

ii) The finite MA(q) process has an ACF that is zero beyond q. For an AR
process, the ACF is infinite in extent and consits of mixture of damped
exponentials and/or damped sine waves.

iii) Parameters of finite MA process are not required to satisfy any

condition for stationarity. However, for invertibility, the roots of the
characteristic equation must lie inside the unit circle.

ARMA modelling is a classic technique which has found a

tremendous number of applications

c Danilo P. Mandic
Digital Signal Processing 46