1. No category

Solutions 3 (STAT673)
(1.8)
(i) Show that the function c(u) = exp(−a|u|) where a > 0 is a positive semi-definite
function.
The Fourier transform is
Z ∞
exp(−a|u|) exp(iuω)du
f (ω) =
−∞
Z ∞
Z 0
exp(−a|u|) exp(iuω)du
exp(−a|u|) exp(iuω)du +
=
−∞
=
a2
0
2a
,
+ ω2
which is clearly positive. Thus c(u) is a valid covariance function.
(ii) Show that the commonly
used exponential spatial covariance defined on R2 ,
p
c(u1 , u2 ) = exp(−a u21 + u22 ), where a > 0, is a positive semi-definite function.
The Fourier transform is
Z ∞Z ∞
q
f (ω1 , ω2 ) =
exp(−a u21 + u22 ) exp(iu1 ω1 + iu2 ω2 )du1 du2
−∞
−∞
Now changing variables from Euclidean to Radial coordinates and using
the Jacobian we have
Z 2π Z ∞
r exp(−ar) exp(ir cos θω1 + ir sin θω2 )drdθ
f (ω1 , ω2 ) =
0
0
Z 2π
1
dθ
=
[i cos θ + i sin θ − a]2
0
2πa
=
(I used Maple)
2
[a + ω12 + ω22 ]3/2
which is non-negative. Thus using the Euclidean distance on c(·) still
gives a valid covariance function.
Xiao also told me that the Hausdorff-Bernstein-Widder Theorem for
completely monotone functions could be used, which would avoid this
calculation.
(2.4) (a) Obtain the stationary solution of the AR(2) process
7
2
Xt = Xt−1 − Xt−2 + εt ,
3
3
where {εt } are iid random variables with mean zero and variance σ 2 .
Does the solution have an MA(∞) representation?
The characteristic function is φ(z) = 1 − (7/3)z + (2/3)z 2 . The roots are
3 and 1/2. Since the roots lie inside and outside the unit circle the
solution is non-casual.
By using partial fractions we have
1
1
6
−
.
=
2
1 − (7/3)z + (2/3)z
5(1 − 2z) 5(1 − z/3)
By using equation (2.5) in the notes this means the solution is
∞
∞
−6 X 1
1X 1
Xt =
ε
−
εt−j .
t+j+1
5 j=0 2j+1
5 j=0 3j
which is in terms of both past and future innovations, thus it is linear
but non-casual.
(b) Obtain the stationary solution of the AR(2) process
√
4× 3
42
Xt =
Xt−1 − 2 Xt−2 + εt ,
5
5
where {εt } are iid random variables with mean zero and variance σ 2 .
Does the solution have an MA(∞) representation?
√
√
φ(z) = 1 − 4 5 3 z + ( 45 )2 , which has roots [(2/5)( 3 ± i)]−1 = (5/4) exp(±iπ/6).
Using partial fractions gives
1
1
(4/5)e−iπ/6
eiπ/6
=
−
(4/5)
.
φ(z)
2(4/5)i sin(π/6)
1 − (4/5)eiπ/6 z (1 − (4/5)e−iπ/6 z)
As the roots lies outside the unit circle the solution is causal we can
use equation (2.4) in the notes to obtain the form
j+1 j+1 !
∞
X
−5i
4 iπ/6
4 −iπ/6
e
Xt =
− e
εt−j
8 sin(π/6) j=0
5
5
∞ j+1
X
5
4
j+1
π εt−j .
=
sin
4 sin(π/6) j=0 5
6
(quite possibly there is something wrong in this calculation!).
(c) Obtain the stationary solution of the AR(2) process
Xt = Xt−1 − 4Xt−2 + εt ,
where {εt } are iid random variables with mean zero and variance σ 2 .
Does the solution have an MA(∞) representation?
The characteristic
function is φ(z) = 1 − z + 4z 2 = (1 − 2e−iθ z)(1 − 2e−θ z),
√
where θ = tan−1 ( 15). Thus we see the roots lie inside the unit circle
and the solution is completely non-causal. Thus the solution is
∞ j
2 X 1
Xt = √
sin (θj) εt+j+1 .
15 j=0 2
(2.5) Construct a causal stationary AR(2) process with pseudo-period 17. Using the R
function arima.sim simulate a realisation from this process (of length 200) and make
a plot of the periodogram. What do you observe about the peak in this plot?
From clas we learnt that to obtain a pseudo-period of period θ we need that
the roots of the characteristic function φ are complex with at argument
2π/17. Thus we let
φ(z) = [1 − zr exp(i2π/17)][1 − zr exp(−i2π/17)] = 1 − 2r cos(
2π
)z + r2 ,
17
where |r| < 1 if the solution is causal (in this case the roots lie outside unit
circle). This corresponds to the model
Xt = 2r cos(
2π
)Xt−1 − r2 Xt−2 + εt .
17
We consider the two cases, r = 4/5 and r = 0.99.
The case r = 4/5 is considered in Figure 1 and the case r = 0.99 is considered
in Figure 2.
Looking at Figure 1a, we observe that in the case r = 4/5, period of 17 is not
so apparent (though is still
This is because we recall the solution
P∞there).
4 j
j)εt . Looking at any one innovation
will have the form Xt = C j=0 ( 5 ) sin( 2π
17
εt and following the coefficient over Xt+j , that is {( 45 )j sin( 2π
j)}, we see that
17
for j = 17 (which is when the new period starts) (4/5)17 = 0.022, hence there
is a lot of damping of the period of the period. This is why the period is
not hugely obvious.
The periodogram and true spectral density is given Figure 1b. There does
seem to be a peak about 2π/17, but it is very irratic, this is because the
periodogram is extremely noisey and it cannot be improved even when the
sample size is large. The smooth version of the periodogram in given below.
Here there is a clear peak at around 2π/17.
In Figure 2 we consider the case r = 0.99 (we are close to the unit root
exp(i2π/17)). We observe that the period now is far more pronounced, this
is because there is less damping in the coefficients of the solution.
5
0
−5
ar2
0
34
68
119
170
221
272
323
374
425
476
Time
0.6
0.2
−0.2
ACF
1.0
Series ar2
0
5
10
15
20
25
Lag
0.20
0.10
0.00
Periodogram
(a) A realisation and the corresponding (empirical) ACF
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
● ●
●●
●
●
●
●
●
●
●
●
●
●● ●
● ●●
●
●
●
●
● ●
●
●
●●
● ●●
●●
●●
●●
● ●● ●
●
●●
●● ●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
● ●● ●
●●●
●●●●●
●●
●●●●●
●●
●
●●
●
●
● ●●●
●●
●
●
●
●
●●
●
●
●
●●●
● ●
●
●
●●●
●
●
●
●
● ●●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
0
1
2
3
4
5
6
frequency
40
20
0
spectrum
60
Autoregressive
0.0
0.1
0.2
0.3
0.4
0.5
frequency
(b) The periodogram and the true spectral density (remember to multiple the x-axis on true spectral
density by 2π to get the x-axis on the periodogram.)
Figure 1: The model Xt = 2 × 0.8 cos( 2π
)Xt−1 − 0.82 Xt−2 + εt
17
(2.6) In the following simulations, use non-Gaussian innovations.
30
10
−30 −10
ar2
0
34
68
119
170
221
272
323
374
425
476
Time
−0.5
ACF
0.5
1.0
Series ar2
0
5
10
15
20
25
Lag
20
●
15
●
●
●
5
10
●
●
0
Periodogram
(a) A realisation and the corresponding (empirical) ACF
●
●
●
●
●
●
●
●
●●
●●
●●
●●
●
●●●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
0
1
2
3
4
5
6
frequency
15000
5000
0
spectrum
Autoregressive
0.0
0.1
0.2
0.3
0.4
0.5
frequency
(b) The periodogram and the true spectral density (remember to multiple the x-axis on true spectral
density by 2π to get the x-axis on the periodogram.)
Figure 2: The model Xt = 2 × 0.99 cos( 2π
)Xt−1 − 0.992 Xt−2 + εt
17
(i) Simulate an AR(4) process with characteristic function
2π
2π
2π
2π
φ(z) = 1 − 0.8 exp(i )z 1 − 0.8 exp(−i )z 1 − 1.5 exp(i )z 1 − 1.5 exp(−i )z .
13
13
5
5
(ii) Simulate an AR(4) process with characteristic function
2π
2π
2
2π
2
2π
φ(z) = 1 − 0.8 exp(i )z 1 − 0.8 exp(−i )z 1 − exp(i )z 1 − exp(−i )z .
13
13
3
5
3
5
Do you observe any differences between these realisations?
I used Furong’s method. I basically simulated an AR(2) process
Xt = 2 × (4/5) × cos(2π/13)Xt−1 − (4/5)2 Xt−2 + εt .
Then to simulate the causal process I used
Yt = 2 × (2/3) × cos(2π/5)Yt−1 − (2/3)2 Yt−2 + Xt
and to simulate the noncausal process I used
Zt−2 = (2/3)2 (−Zt + 2 × (3/2) × cos(2π/5)Zt−1 + Xt ) .
The code is given at the end of the solutions.
The plot of the causal and noncausal realisations are given in Figure 3
The corresponding ACF plot is given in Figure 4. We observe that they
are identical, which is what we would expect. However, when I calculate
their skewness and kurtosis we observe something very interesting. The
skewness of the causal process is (from 300-600) is −0.13 and the non-causal
it is −0.189 (not much difference there). However, the kurtosis for the causal
time series is 1.47 whereas the kurtosis for the non-causal time series is 0.48
(which is a considerable difference). Here we are seeing differences in the
higher order cumulants.
15
5
−5
−15
Y[c(300:600)]
300
350
400
450
500
550
600
500
550
600
2 4 6
−2
−6
Z[c(300:600)]
c(300:600)
300
350
400
450
c(300:600)
Figure 3: The top plot is the causal realisation and the lower plot is the non-causal realisation.
0.2
−0.2
ACF
0.6
1.0
Series Y[c(300:600)]
0
5
10
15
20
Lag
0.2
−0.2
ACF
0.6
1.0
Series Z[c(300:600)]
0
5
10
15
20
Lag
Figure 4: The top plot is the acf of causal realisation and the lower plot is the acf non-causal
realisation (using the observations from 300-600).
X = c (1 ,1000)
Y = c (1 ,1000)
Z = c (1 ,1000)
# S i m u l a t e from a t−d i s t r i b u t i o n with 3 d f
rt3 = rt (1000 ,3)
# I n i t i a l values
X[ 1 ] = rt3 [ 1 ]
X[ 2 ] = rt3 [ 2 ]
Y[ 1 ] = rt3 [ 1 ]
Y[ 2 ] = rt3 [ 2 ]
# Causal
for ( i in c (3:1000)){
X[ i ] = 2 ∗ 0 . 8 ∗ c o s (2∗ p i /13)∗X[ i −1] − ( 0 . 8 ∗ ∗ 2 ) ∗X[ i −2] + r t 3 [ i ]
Y[ i ] = 2 ∗ ( 2 / 3 ) ∗ c o s (2∗ p i /5)∗Y[ i −1] − ( ( 2 / 3 ) ∗ ∗ 2 ) ∗Y[ i −2] + X[ i ]
}
# Noncausal ( s e t t i n g i n i t i a l v a l u e s )
Z [ 1 0 0 0 ] = X[ 1 0 0 0 ]
Z [ 9 9 9 ] = X[ 9 9 9 ]
for ( i in c (1000:3)){
Z [ i −2] = ((2/3) ∗∗2)∗( −Z [ i ] + 2 ∗ ( 3 / 2 ) ∗ c o s (2∗ p i /5)∗Z [ i −1] + X[ i ] )
}
# To measure k u r t o s i s and ske wne ss
k u r t o s i s 1 <− f u n c t i o n ( x ) {
m4 <− mean ( ( x−mean ( x ) ) ˆ 4 )
k u r t <− m4/ ( sd ( x)ˆ4) −3
kurt
}
s k e w n e s s 1 <− f u n c t i o n ( x ) {
m3 <− mean ( ( x−mean ( x ) ) ˆ 3 )
skew <− m3/ ( sd ( x ) ˆ 3 )
skew
}
# P l o t s ( u s i n g burn−i n from e i t h e r s i d e s )
par ( mfrow=c ( 2 , 1 ) )
p l o t ( c ( 3 0 0 : 6 0 0 ) ,Y[ c ( 3 0 0 : 6 0 0 ) ] , type=” l ” )
p l o t ( c ( 3 0 0 : 6 0 0 ) , Z [ c ( 3 0 0 : 6 0 0 ) ] , type=” l ” )