Genshiro Kitagawa. Introduction to Time Series Modeling (Введение в моделирование временных рядов)
Подождите немного. Документ загружается.
Append ix D
Algorithm for the Monte Carlo Filter
This a ppendix pr esents a derivation of the Monte Carlo filter algorithm.
The readers are referred to Kitagawa (1996) and Doucet et al. (2001) for
details.
D.1 One-Step- A head Prediction
Assume that m particles {f
(1)
n−1
,.. ., f
(m)
n−1
} that can be considered as m
indepen dent realizations from the co nditional distribution p(x
n−1
|Y
n−1
)
of the state x
n−1
, and m particles {v
(1)
n
,.. ., v
(m)
n
} fro m the system noise
density p(v) are given. Namely, assume that
f
(i)
n−1
∼ p(x
n−1
|Y
n−1
), v
(i)
n
∼ q(v). (D.1)
Then, it can be shown that the one-step-ahead predic tive distribution
of x
n
is given by
p(x
n
|Y
n−1
)
=
Z Z
p(x
n
,x
n−1
,v
n
|Y
n−1
)dx
n−1
dv
n
(D.2)
=
Z Z
p(x
n
|x
n−1
,v
n
,Y
n−1
)p(v
n
|x
n−1
,Y
n−1
)p(x
n−1
|Y
n−1
)dx
n−1
dv
n
.
Since the system noise v
n
is independent of past states and ob-
servations, the conditio nal distribution of the system noise satisfies
p(v
n
|x
n−1
,Y
n−1
) = p(v
n
). On the other hand, since x
n
depends only on
x
n−1
and v
n
, p(x
n
|x
n−1
,v
n
,Y
n−1
) = p(x
n
|x
n−1
,v
n
) =
δ
(x
n
−F(x
n−1
,v
n
)),
and we have
p(x
n
|Y
n−1
) =
Z Z
δ
(x
n
−F(x
n−1
,v
n
))p(v
n
)p(x
n−1
|Y
n−1
)dx
n−1
dv
n
.
(D.3)
Therefore, when realizations {v
( j)
n
}of p(v
n
) and {f
( j)
n−1
}of p(x
n−1
|Y
n−1
)
are given, realizations {p
( j)
n
} of p(x
n
|Y
n−1
) are obtained by
p
( j)
n
= F( f
( j)
n−1
,v
( j)
n
). (D .4)
259
260 ALGORITHM FOR THE MONTE CARLO FILTER
D.2 Filter
When m independen t realiza tions p
(1)
n
,.. ., p
(m)
n
of the distribution
p(x
n
|Y
n−1
) are given, it is equivalent to approximate the distribution
p(x
n
|Y
n−1
) by the emp irical distribution function
P
n
(x) =
1
m
m
∑
i=1
I(x, p
(i)
n
), (D.5)
where I(x,a) = 0 for x < a and I(x,a) = 1 otherwise. This means that
the predictive distribution p (x
n
|Y
n−1
) is approximated by the probability
function
Pr(x
n
= p
( j)
n
|Y
n−1
) =
1
m
, for j = 1,.. .,m. (D.6)
Then, given the observation y
n
, the poster ior distribution of x
n
is ob-
tained by
Pr(x
n
= p
( j)
n
|Y
n
) = Pr(x
n
= p
( j)
n
|Y
n−1
,y
n
)
= lim
∆y→0
Pr(x
n
= p
( j)
n
,y
n
≤ y ≤ y
n
+ ∆y|Y
n−1
)
Pr(y
n
≤ y ≤ y
n
+ ∆y|Y
n−1
)
=
p(y
n
|p
( j)
n
)Pr(x
n
= p
( j)
n
|Y
n−1
)
m
∑
i=1
p(y
n
|p
(i)
n
)Pr(x
n
= p
(i)
n
|Y
n−1
)
=
α
( j)
n
·
1
m
m
∑
i=1
α
(i)
n
·
1
m
=
α
( j)
n
m
∑
i=1
α
(i)
n
. (D.7)
The cumula tive distribution f unction
1
m
∑
i=1
α
(i)
n
m
∑
i=1
α
(i)
n
I(x, p
(i)
n
) (D.8)
correspo nding to this Pr(x
n
= p
( j)
n
|Y
n
) has jumps with jump sizes pro-
portion al to
α
(1)
n
,.. .,
α
(m)
n
given by the right- hand side of the equ ation
(15.9) only at the m p oints p
(1)
n
,.. ., p
(m)
n
.
Although the approximation of the distribution of the filter was ob-
tained using this expression (D.8), it is convenient to re-approximate it
SMOOTHING 261
by m par ticles f
(1)
n
,.. ., f
(m)
n
with equal weights to perform the compu-
tation of the prediction (D.1) in the next time step. This corresponds to
representin g the distribution of (D.8) by the emp irical distribution func-
tion
1
m
m
∑
i=1
I(x, f
(i)
n
). (D.9)
The m rea lizations {f
(1)
n
,.. ., f
(m)
n
} can be obtained by re-sampling
{p
(1)
n
,.. ., p
(m)
n
} with proba bilities
Pr( f
( j)
n
= p
(i)
n
|Y
n
) =
α
(i)
n
α
(1)
n
+ ···+
α
(m)
n
, j = 1, ... ,m. (D.10)
D.3 Smoothing
For smoothing, assume that Pr(x
1
= s
( j)
1|n−1
,.. ., x
n−1
= s
( j)
n−1|n−1
|Y
n−1
) =
1/m and v
( j)
n
∼ q(v) and define (p
( j)
1|n−1
,.. ., p
( j)
n|n−1
) as follows:
p
( j)
i|n−1
=
s
( j)
i|n−1
, for i = 1,. .., n −1
F(s
( j)
n−1|n−1
,v
( j)
n
), for i = n.
(D.11)
Then, (p
( j)
1|n−1
,.. ., p
( j)
n|n−1
) can be considered as a realization from the
joint distribution of (x
1
,.. ., x
n
) when the observation Y
n−1
is given.
Next, given the observation y
n
, the distribution Pr(x
1
≤ p
( j)
1|n−1
,.. ., x
n
≤
p
( j)
n|n−1
|Y
n−1
) can be updated as fo llows:
Pr(x
1
= p
( j)
1|n−1
,.. ., x
n
= p
( j)
n|n−1
|Y
n
)
= Pr(x
1
= p
( j)
1|n−1
,.. ., x
n
= p
( j)
n|n−1
|Y
n−1
,y
n
)
=
p(y
n
|x
1
= p
( j)
1|n−1
,.. ., x
n
= p
( j)
n|n−1
)
p(y
n
|Y
n−1
)
× Pr(x
1
= p
( j)
1|n−1
,.. ., x
n
= p
( j)
n|n−1
|Y
n−1
)
=
p(y
n
|p
( j)
n|n−1
)Pr(x
1
= p
( j)
1|n−1
,.. .,x
n
= p
( j)
n|n−1
|Y
n−1
)
p(y
n
|Y
n−1
)
.(D.12)
Since p
( j)
n|n−1
is the same as the particle p
( j)
n
of the filter algorithm
262 ALGORITHM FOR THE MONTE CARLO FILTER
(D.5), the smoothing distribution p(x
1
,.. ., x
n
|Y
n
) can be obtained by
re-sampling m n-dimensional vectors (p
( j)
1|n−1
,.. ., p
( j)
n|n−1
)
T
, j = 1,.. .,m
with the same weights as the filter.
263
Answers to the Problems
Chapter 1
1. It is important to select the sampling in te rval appropriately. If it is too
wide, it is not possible to capture the feature s of the continuous time
series. On the oth e r hand, if it is too narrow, many parameters in the
modeling may be used up o n nonessentials, thus preventing effective
representation of the time series. A rough r ule of thumb is to set the
sampling interval as 1/5 to 1/10 of the dominant period.
2. (An example.) Time series obtained by r e cording the sales amount of
instant coffee at a store. Since values are non-negative, the distribu-
tion of the time series is asymmetric. Also this series may suddenly
increase due to discounting o r due to promo tion activities.
3.(1) By solving (1.1) with respe c t to y, we obtain y = e
z
/(1 + e
z
).
(2) z = log{(y −a)/(b −y)}, and the inverse transformation is given
by y = (a + be
z
)/(1 + e
z
).
4. If the time series contains ob servation noise and is expressed as
y
n
= a + bn +
ε
n
, the difference of y
n
becomes ∆y
n
= b +
ε
n
−
ε
n−1
.
Therefore, although th e trend component is removed, the noise com-
ponen ts become more complex.
5. The corrected value of this y ear is affected by the change of the trend
in the previous year. For example, if the slope of the tr end incr e ased in
the middle of the previous year, the annual rate looks as if it decreased
from the middle of this year.
6.(1) Assume that y
n
= T
n
+ w
n
, T
n
= a + bn, and w
n
∼ N(0,
σ
2
) is a
white noise, then
ˆ
T
n
=
1
3
(y
n−1
+ y
n
+ y
n+1
)
=
1
3
(T
n−1
+ T
n
+ T
n+1
) +
1
3
(w
n−1
+ w
n
+ w
n+1
).
Here, the first ter m on the right-hand side is T
n
. On the other
hand, the mean of the second term is 0 and fr om E(w
n−1
+ w
n
+
w
n+1
)
2
= 3
σ
2
, the variance becomes
σ
2
/3.
(2) By setting the number of terms large, we c an get smoother esti-
mates. On the o ther hand, this makes it difficult to detect sudden
structural c hanges and estima te s may become sensitive to outlying
observations. The m oving median has th e opposite properties.
264 ANSWERS TO THE PROBLEMS
Chapter 2
1. A Gaussian distribution is com pletely specified by its mean and vari-
ance. Theref ore, if the mean and the variance are time-invariant,
this also means that the distribution is tim e -invariant an d becomes
strongly stationary.
2. Co nsider a time series with a standard Cauchy distribution. Then ob-
viously, it is stron gly stationary. However, since the Cauchy distribu-
tion does not have a mean or variance, th e series cannot be weakly
stationary.
3. C
k
= E(y
n
−
µ
)(y
n+k
−
µ
) = E(y
n+k
−
µ
)(y
n
−
µ
) = E(y
n
−
µ
)(y
n−k
−
µ
) = C
−k
. The third equality holds only for a stationary
time series y
n
.
4. C
0
= Ey
2
n
= Ev
2
n
− 2cEv
n
v
n−1
+ c
2
Ev
2
n−1
= 1 + c
2
, C
1
= E(v
n
−
cv
n−1
)(v
n−1
−cv
n−2
) = −c, C
k
= 0 (k ≥ 2).
5. For arbitrary
α
= (
α
1
,.. .,
α
k
) with not all co mponents ze ro,
α
C
α
T
=
k
∑
i=1
k
∑
j=1
α
i
α
j
C
i−j
=
k
∑
i=1
k
∑
j=1
α
i
α
j
E(y
n−i
y
n−j
)
= E
k
∑
i=1
α
i
y
n−i
2
≥ 0
6.(1) By taking the expectation of both sides,
E[
ˆ
C
k
(i, j )] =
1
N
N
∑
n=k+1
C
k
(i, j ) =
N −k
N
C
k
(i, j )
(2) Since the sample autocovariance function defined by (2.21) is pos-
itive semi-definite, it has the significant advantage that the esti-
mated AR models always become stationary. If N is replaced by
N −k, we can get an unbiased estimate of C
k
(i − j). Instead how-
ever, we lose the advantage above of stationarity.
7.(1) When the sample size is n,
ˆ
C
0
∼ N(0, 2
σ
4
/n),
ˆ
C
k
∼ N(0,
σ
4
/n),
ˆ
R
k
∼ N(0,n
−1
).
(2) For all k, chec k if |
ˆ
R
k
|<
α
√
n
holds. For example.
α
can be set equal
to 2.
ANSWERS TO T H E PROBLEMS 265
Chapter 3
1. From the expression e
−2
π
ik f
= cos(2
π
k f ) −isin(2
π
k f ), we have
p( f ) =
∞
∑
k=−∞
C
k
cos(2
π
k f ) −i
∞
∑
k=−∞
C
k
sin(2
π
k f ) . Here , since C
k
is
an even function and sin(2
π
k f ) is an odd function, C
k
sin(2
π
k f ) be-
comes an odd function and the seco nd term on the right-hand side of
the above eq uation becomes 0.
2. p( f ) = 1 +C
2
−2Ccos(2
π
f )
3.
p( f ) =
∞
∑
k=−∞
C
k
e
−2
π
ik f
=
σ
2
1 −a
2
∞
∑
k=−∞
a
|k|
e
−2
π
ik f
=
σ
2
1 −a
2
1
1 −ae
−2
π
i f
+
1
1 −ae
2
π
i f
−1
=
σ
2
1 −a
2
1 −a
2
|1 −ae
−2
π
i f
|
2
=
σ
2
|1 −ae
−2
π
i f
|
2
4.
1
N
N
∑
n=1
y
n
exp(−2
π
i(n −1) j/N)
2
=
1
N
N
∑
n=1
N
∑
m=1
y
n
y
m
exp(−2
π
i(n −1) j/N)exp(2
π
i(m −1) j/N)
=
1
N
(
N
∑
n=1
y
2
n
+ 2
N−1
∑
k=1
N
∑
n=k+1
y
n
y
n−k
exp(−2
π
ik j/N)
)
=
ˆ
C
0
+ 2
N−1
∑
k=1
ˆ
C
k
exp(−2
π
ik j/N)
5. By taking the expec ta tion of both sides o f (3.9),
E[ ˆp( f )] =
N−1
∑
k=1−N
E[
ˆ
C
k
]e
−2
π
ik f
=
N−1
∑
k=1−N
N −|k |
N
C
k
e
−2
π
ik f
=
N−1
∑
k=1−N
C
k
e
−2
π
ik f
−
1
N
N
∑
k=1−N
|k|C
k
e
−2
π
ik f
266 ANSWERS TO THE PROBLEMS
Chapter 4
1. The log-likelihood is given by
ℓ =
n
∑
i=1
log f (m
i
|
λ
) =
n
∑
i=1
{−
λ
+ m
i
log
λ
−logm
i
!}.
Therefore, from
∂
ℓ
∂λ
=
n
∑
i=1
n
−1 +
m
i
λ
o
= −n +
1
λ
n
∑
i=1
m
i
= 0
we have
ˆ
λ
= n
−1
n
∑
i=1
m
i
.
2.(1)
ˆ
µ
0
=
1
n + m
n
∑
i=1
x
i
+
m
∑
i=1
y
i
!
,
ˆ
µ
1
=
1
n
n
∑
i=1
x
i
,
ˆ
µ
2
=
1
m
m
∑
i=1
y
i
,
ˆ
σ
2
0
=
1
n + m
n
∑
i=1
(x
i
−
ˆ
µ
0
)
2
+
m
∑
i=1
(y
i
−
ˆ
µ
0
)
2
!
,
ˆ
τ
2
0
=
1
n + m
n
∑
i=1
(x
i
−
ˆ
µ
1
)
2
+
m
∑
i=1
(y
i
−
ˆ
µ
2
)
2
!
.
Then by comp aring
AIC
0
= (n + m)(log2
π
ˆ
σ
2
0
+ 1) + 2 ×2
AIC
1
= (n + m)(log2
π
ˆ
τ
2
0
+ 1) + 2 ×3
and if AIC
0
< AIC
1
, then this means that we can consider the two
means to be identical.
(2)
ˆ
σ
2
1
=
1
n
n
∑
i=1
(x
i
−
ˆ
µ
0
)
2
,
ˆ
σ
2
2
=
1
m
m
∑
i=1
(y
i
−
ˆ
µ
0
)
2
, AIC
2
= n(log 2
π
ˆ
σ
2
1
+
1) + m(log 2
π
ˆ
σ
2
2
+ 1) + 2 × 3. If AIC
0
< AIC
2
, then the two
variances can be considered identical. Here, from
∑
∞
k=−∞
|C
k
| <
∞, for arbitrary
ε
> 0, there exists an integer m tha t satisfies
∑
∞
k=m+1
|C
k
| <
ε
/2. Then, we have
1
N
∞
∑
k=−∞
|k|C
k
e
−2
π
ik f
≤
2
N
∞
∑
k=1
k|C
k
| ≤
2
N
m
∑
k=1
kC
0
+ 2
∞
∑
k=m+1
k
N
|C
k
|
≤
2
N
m(m + 1)
2
C
0
+
ε
ANSWERS TO T H E PROBLEMS 267
Since by taking N sufficiently large, the first term on the right-hand
side of the equation can be made arbitrarily small, it follows that
E [ ˆp( f )] → p( f ).
3. Putting
ˆ
µ
=
1
n
n
∑
i=1
y
i
,
ˆ
σ
2
=
1
n
n
∑
i=1
(y
i
−
ˆ
µ
)
2
, the value of the AIC for
the Gaussian distribution model is obtained by AIC
0
= n log 2
π
ˆ
σ
2
+
n + 2 ×2. On the o ther hand, the AIC value for the Cauchy model is
obtained by AIC
1
= n log
ˆ
τ
2
−2n log
π
−
∑
n
i=1
log{(y
i
−
ˆ
µ
)
2
+
ˆ
τ
2
}+
2 ×2, where
ˆ
µ
and
ˆ
τ
2
are obtained by m aximizing (4.36). T hen, if
AIC
0
< AIC
1
, the Gaussian model is considered to be better.
4.(1) For the binomial distribution model f (m|n, p) =
n
C
m
p
m
(1 −
p)
n−m
, the log-likelihood is defined by ℓ(p) = log
n
C
m
+ m log p +
(n −m)log(1 − p). Therefore, from
dℓ
d p
= 0,
m
p
−
n −m
1 − p
= 0, we
obtain ˆp = m/n.
(2) The AIC value for the model (1) is obtained by AIC
1
=
−2 lo g
n
C
m
− 2mlog ˆp − 2(n − m)log(1 − ˆp ) + 2. On the o ther
hand, if we assume equal pro babilities, taking p = 1/2, we ob-
tain AIC
0
= −2log
n
C
m
+ 2nlog2. Here, if AIC
0
< AIC
1
, we can
consider that the equal probability m odel is better.
5. From
d
d
θ
log f (y|
θ
) = f (y|
θ
)
−1
d f (y|
θ
)
d
θ
, we have
d
2
d
θ
2
log f (y|
θ
) =
d
d
θ
d
d
θ
log f (y|
θ
)
=
d
d
θ
f (y|
θ
)
−1
d f
d
θ
= −f (y|
θ
)
−2
d f
d
θ
2
+ f (y|
θ
)
−1
d
2
f
d
θ
2
= −
d
d
θ
log f (y|
θ
)
2
+ f (y|
θ
)
−1
d
2
f
d
θ
2
.
Therefore,
J = −E
d
2
d
θ
2
log f (y|
θ
)
= E
d
d
θ
log f (y|
θ
)
2
−E
f (y|
θ
)
−1
d
2
f
d
θ
2
268 ANSWERS TO THE PROBLEMS
= I −
Z
∞
−∞
d
2
f
d
θ
2
dy = I −
d
2
d
θ
2
Z
∞
−∞
f (y|
θ
)dy = I
6. From z = k(y) = logy, f (z) = (2
πσ
2
)
−1
exp{−(z −
µ
)
2
/2
σ
2
},
and we have
dk
dy
= y
−1
. Ther efore,
g(y) =
dk
dy
f (k(y)) =
1
y
1
√
2
πσ
2
e
−
(log y−
µ
)
2
2
σ
2
.
Chapter 5
1. Let
ˆ
µ
denote the maximum likelihood estimate of
µ
. Since the num-
ber of parameters is 1, the AIC is obtained by
AIC = n log 2
πσ
2
−
1
σ
2
n
∑
i=1
(y
i
−
ˆ
µ
)
2
+ 2.
2.(1) From S =
∑
N
n=1
(y
n
−ax
2
n
−bx
n
)
2
,
dS
da
= −2
N
∑
n=1
x
2
n
(y
n
−ax
2
n
−bx
n
)
= −2
N
∑
n=1
x
2
n
y
n
−a
N
∑
n=1
x
2
n
−b
N
∑
n=1
x
3
n
= 0
dS
db
= −2
N
∑
n=1
x
n
(y
n
−ax
2
n
−bx
n
)
= −2
N
∑
n=1
x
n
y
n
−a
N
∑
n=1
x
3
n
−b
N
∑
n=1
x
2
n
= 0.
Solving these, we have
ˆa =
∑
x
2
n
y
n
∑
x
2
n
−
∑
x
n
y
n
∑
x
3
n
∑
x
4
n
∑
x
2
n
−(
∑
x
3
n
)
2
ˆ
b =
∑
x
n
y
n
∑
x
4
n
−
∑
x
2
n
y
n
∑
x
3
n
∑
x
4
n
∑
x
2
n
−(
∑
x
3
n
)
2
.
(2) y
n
= ax
n
(x
n
−b) +
ε
n
. By minimiz ing S =
N
∑
n=1
(y
n
−ax
2
n
−abx
n
)
2
with respect to the two parameters a and b by a num erical opti-
mization proced ure, we can ob ta in the least squares estimates.