SDS-1

$5REXVWILOWHULQJ0HWKRG for Noisy
Non-6tationary 0XOWLYDULDWHTime
SeriesZLWK(FRQRPHWULF$SSOLFDWLRQV
1DRWR.XQLWRPR DQGSeisho Sato
2FWREHU 2020

Statistics & Data Science Series back numbers:
http://www.mims.meiji.ac.jp/publications/datascience.html

A Robust-ﬁltering Method for Noisy
Non-Stationary Multivariate Time Series with
Econometric Applications ∗Naoto Kunitomo †andSeisho Sato ‡October 16, 2020
(A Revised Version)

Abstract
We investigate a new ﬁltering method to estimate the hidden states of randomvariables for multiple non-stationary time series data. This helps in analyzing
small sample non-stationary macro-economic time series in particular and it isbased on the frequency domain application of the separating information maximum likelihood (SIML) method, developed by Kunitomo, Sato and Kurisu(2018), and Kunitomo, Awaya and Kurisu (2019), and Nishimura, Sato and
Takahashi (2019). We solve the ﬁltering problem of hidden random variablesof trend-cycle, seasonal and measurement-errors components, and propose amethod to handle macro-economic time series. We develop the asymptotic
theory based on the frequency domain analysis for non-stationary time series.
We illustrate applications, including some properties of the method of Müller
and Watson (2018), and analyses of some macro-economic data in Japan.

∗

This is a revised version of the paper (SDS-9, MIMS=RSP Statistics and Data Science Series,
Meiji University) presented at the conference ”Theory and Applications of Economic and Government Statistics 2019” held at University of Tokyo on January 31, 2019. We thank Akihiko
Takahashi for giving a valuable suggestion on the SIML ﬁltering method at the initial stage of thiswork.
†
The corresponding author, School of Political Science and Economics, Meiji University, Kanda-Surugadai 1-1 (Sarugakucho 3rd Building C-106), Chiyoda-ku, Tokyo 101-8301,
naoto.kunitomo@gmail.com‡
Graduate School of Economics, University of Tokyo, Bunkyo-ku,Hongo7-3-1, Tokyo113-0033,
seisho@e.u-tokyo.ac.jp

1

Key Words
Noisy Non-stationary time series, Errors-variables models, Measurement Error, trendcycle, seasonality and noise, Robust-ﬁltering, SIML, Real-Fourier-Inversion, MüllerWatson method, Macro-consumption in Japan

2

1

Introduction

There exists vast research on the use of statistical time series analysis for macroeconomic time series. One important feature of macroeconomic time series, which
is diﬀerent from standard statistical time series analysis, is that the observed timeseries is an apparent mixture of non-stationary and stationary components. The
second feature is that the measurement errors in economic time series play important roles, because macro-economic data are usually constructed from various
sources, including sample surveys in major oﬃcial statistics, while the statisticaltime series analysis often ignored measurement errors. Third, the sample size ofmacro-economic data is rather small, and we have about 120 time series observations for each series for quarterly data over 30 years. The quarterly GDP series,
which is the most important data in the Japanese macro-economy, for instance, isregularly estimated and released by the cabinet oﬃce of Japan 1 . Fourth, to publish
the seasonally adjusted data, the oﬃcial agencies usually apply the X-12-ARIMA
program, which uses the univariate reg-ARIMA model to remove the seasonality as
the standard ﬁltering procedure. As the sample size is small, it is important to use
an appropriate statistical procedure to extract information on trend-cycle, seasonality and noise (or measurement error) components systematically from multiple timeseries data.
In this study, we investigate a new ﬁltering procedure to estimate the hidden
states of trend-cycle, which are non-stationary, and to handle multiple time series
data, including small sample time series. Kunitomo and Sato (2017), Kunitomo,
Sato and Kurisu (2018), and Kunitomo, Awaya and Kurisu (2019) have developed
the separating information maximum likelihood (SIML) method for estimating thenon-stationary errors-in-variables models. They have discussed the asymptotic andﬁnite sample properties of the estimation of unknown parameters in the statisticalmodels. We utilize their results to solve the ﬁltering problem of hidden randomvariables and show that they lead to new a way to handle macro-economic timeseries.
Related literature on the non-stationary economic time series analysis are Engle and Granger (1987) and Johansen (1995), which examined multivariate nonstationary and stationary time series and developed the notion of co-integration
without measurement errors. Our problem is related to their work, but it has different aspects, and our focus is on the non-stationary trend-cycle, seasonality andmeasurement error in the non-stationary errors-in-variable models and their frequency domain analysis. Some related econometric studies on time series in the1

Currently, GDP time series data in Japan are constructed and seasonally adjusted since
1994Q1, although some may think the historical data before 1994 was constructed and measured
exactly in the same way. The measuring procedures of oﬃcial GDP series has been changed several times including the base-year changes. See https://www.esri.cao.go.jp/index-e.html for thedetailed explanation of how GDP in Japan is constructed.

3

frequency domain are Baxter and King (1999), Christiano and Fitzgerald (2003)
and Müller and Watson (2018). See Yamada (2019) for a survey of related studies,including the well-known Hodrick-Prescot (HP) ﬁlter in econometrics.
In statistical multivariate analysis, some studies on the errors-in-variables models
are Anderson (1984, 2003) and Fuller (1987); however, they considered multivariate
cases of independent observations, and the underlying situation is diﬀerent fromours.
For statistical ﬁltering methods, Kitagawa (2010) discussed the standard statistical methods already known, including the Kalman-ﬁltering and particle-ﬁltering
methods. Although many studies have examined statistical ﬁltering theories, wemust exercise caution in analyzing non-stationary multivariate economic time series. See Granger and Hatanaka (1964), Brillinger and Hatanaka (1969) on earlystudies, and Harvey and Trimbur (2008) on the relationship between HP ﬁlter andother methods, for instance. Here we should mention two issues. First, the existing
methods often depend on the underlying distributions such as the Gaussian distributions for the Kalman-ﬁltering, and second, they often depend on the dimensionof state variables. There may be some diﬃculty in extending the existing methods
to high-dimension cases, even when the dimension is about 10. On the other hand,
we expect that our method has robustness properties when we handle small sample
economic times series with non-stationary trend-cycle, and stationary seasonality
and measurement errors because our method does not depend on the speciﬁc distribution as well as the dimension of the underlying random variables. See Kunitomo,
Awaya and Kurisu (2017) for a comparison of small sample properties of the ML(maximum likelihood) and SIML methods for the non-stationary errors-in-variables
models, and Nishimura, Sato and Takahashi (2019) for an application of ﬁnancial data smoothing. The most important feature of the present procedure is thatit may be applicable to small sample time series data with non-stationary trendcycle, seasonal and noise components and it has a statistical foundation based onthe (real-valued) spectral decomposition of stochastic processes by a (real-valued)
Fourier transformation, as we shall explain in Sections 4 and 5.
In Section 2, we give some macro-economic data, which have motivated thispresent study. In Section 3, we deﬁne the non-stationary errors-in-variables models
and the SIML method. In Section 4, we introduce the SIML ﬁltering method. In
Section 5, we give the statistical foundation of the method and in Section 6, wediscuss the problem of choosing the number of orthogonal processes, and give somenumerical examples based on simulation and data. Section 7 contains some applications including an interpretation on the Müller-Watson method in econometrics
and gives two empirical applications of macro-consumption in Japan. Concluding
remarks are given in Section 8. The Appendix contains some mathematical derivations of our results and ﬁgures.

4

2

Two Illustrative Examples

In the ﬁrst illustrative example, we plot the graph of two macro-economic time series
in Japan: quarterly (real) consumption and quarterly (real) GDP (1994Q1-2018Q2)
as Figure 1 2 . It looks like a simple example of linear regression in undergraduate
textbooks. However, if we draw the time series sequences of these (original andseasonally unadjusted oﬃcial) macro-data estimated by the Cabinet oﬃce of theJapanese Government as Figure 2, we ﬁnd things to be a little more complex thanin Figure 1. There are clear trend-cycle components, seasonal ﬂuctuations, andnoise components in two time series data. Although many economists usually use
the seasonally adjusted (published) data, which were constructed by using the X-12ARIMA program in diﬀerent ministries within the government, the eﬀects of ﬁlteringin the program are often unknown. The X-12-ARIMA program uses the univariate
reg-ARIMA model, which is a mixture of univariate seasonal ARIMA and linear
regression, and it decomposes univariate time series into the trend-cycle, seasonality
and noise components as the standard procedure 3 . In contrast, the DECOMP
program, which is explained by Kitagawa (2010), is a possible choice particularly inJapan, uses the univariate AR model and the Kalman ﬁltering technique with AIC,
which is based on the Gaussian likelihood. When each time series is handled using
diﬀerent ﬁltering procedures (i.e., diﬀerent ARIMA models or reg-ARIMA models
for instance), it may cause a fundamental problem in their interpretation when thefocus is on the relationships among diﬀerent non-stationary time series.
Figure 3 gives three diﬀerent macro-consumption data (2002 January - 2016 December), which are observed as monthly time series and widely used by economists inJapan to judge the current macro-business condition. The ﬁrst series is Kakei-Chosa(the data from monthly consumer-survey collected by Statistics Bureau, Ministryof Internal Aﬀairs and Communications), the second is Shougyo-Doutai-Statistics(the data from monthly retail constructed by Ministry of Economy, Trade and Industry (METI)). and the third is Dai-Sanji-Sangyo-Statistics (the index data oncommerce constructed by METI). We note that the data construction processesof these series based on sample surveys are complex in diﬀerent ministries within
the government, and each data reﬂects diﬀerent aspects of macro-consumption. Although they show similar movements, we observe some diﬀerences in trend-cycle,
seasonality and noises. Then it may be desirable to unify the monthly consumption
series because we want to judge the business condition each month by just observing
these data to evaluate the state of the Japanese macro-economy and forming macroeconomic policy. Many economists in both governments and private sectors usually
use seasonally adjusted data, which were constructed from the quarterly or monthly(original) time series via the univariate X-12-ARIMA seasonal adjustment program.
2

In Japan, both the original quarterly series and the seasonally adjusted series of GDP and its
major components are regularly published. It diﬀers from macro-data in the U.S. in some aspect.
3
See Findley et.al. (1998) for the detail of X-12-ARIMA program.

5

11.30
11.25
11.20
11.15

RCONSUMPTION(log)

11.10
11.05

11.60

11.65

11.70

11.75

11.80

RGDP(log)

125000110000RGDP

Figure 1: GDP vs. Consumption. (Data are the Quarterly real GDP and realconsumption between 1994Q1-2014Q3, which were published in 2015 by the Economic Social
Research Institute (ESRI), Cabinet Oﬃce, Japan.)

1995

2000

2005

2010

2015

2010

2015

75000
65000

RCON

Time

1995

2000

2005
Time

Figure 2: GDP and Consumption. (Data are the Quarterly real GDP and realconsumption between 1994Q1-2014Q3, which were published in 2015 by the Economic Social
Research Institute (ESRI), Cabinet Oﬃce, Japan.)6

Figure 3: Monthly Consumption Series. (Data are the monthly consumptions between
2002M1-2016M12, which were published in 2017 by the Statistics Bureau and Ministry of Economy,
Trade and Industry (METI), Japan.)
It is important to construct the monthly consumption index, which is consistent
with the published quarterly macro-consumption data, which usually reported with
substantial time lags. It was one of motivations to develop our ﬁltering theory.
Some econometricians use the multivariate (parametric) time series models such
as VAR (vector autoregressive process) for analyzing macro-economic data and investigating the relationships among them. They may use the seasonally adjusted(oﬃcial) data, but we need caution to use such data because most published oﬃcial data are already ﬁltered by the X-12-ARIMA program. When the dimension
is more than 2, some diﬃculties handling trend-cycle, seasonal, and measurement
errors at the same time typically arise. A need to handle macro-economic data insimple non-parametric way motivated us to develop the multivariate non-stationary
errors-in-variables models and the ﬁltering method for the hidden state variables
with measurement errors in Sections 3 and 4.
7

3

Non-Stationary Errors-in-Variables Model and
SIML

3.1

A Simple Non-Stationary Errors-in-Variables Model andthe SIML Method

In this sub-section we ﬁrst introduce a simple non-stationary errors-in-variables
model with trend and noise components, and explain the SIML method. Thenin the next subsection we shall investigate the general framework for non-stationary
multivariate time series with trend-cycle, seasonal and measurement-error components.
Let yji be the i−th observation of the j−th time series at i for i = 1, · · · , n; j =1, · · · , p. We set yi = (y1i , · · · , ypi ) be a p × 1 vector and Yn = (yi ) (= (yij )) bean n × p matrix of observations and denote y0 as the initial p × 1 vector and it
is ﬁxed. We consider the simple model that the underlying non-stationary trend
is xi (= (x1i , · · · , xpi ), i = 1, · · · , n) and the noise component is vi = (v1i , · · · , vpi ),which is independent of xi . We write
yi = xi + vi (i = 1, · · · , n) .

(3.1)

When each pair of vectors Δxi and vi are independently, identically, and normally
distributed (i.i.d.) as Np (0, Σx ) and Np (0, Σv ), respectively, we have the observa tions of an n × p matrix Yn = (yi ) and set the np × 1 random vector (y1 , · · · , yn ) .
Given the initial condition y0 (= x0 ), we have




(3.2)
vec(Yn ) ∼ Nn×p 1n · y0 , In ⊗ Σv + Cn Cn ⊗ Σx ,
where 1n = (1, · · · , 1) and

⎛

1 0
⎜ 1 1
⎜
Cn = ⎜⎜ 1 1⎝ 1 ···
1 ···

(3.3)

···
0
1
1
1

0
···
···
1
1

⎞
0
0 ⎟
⎟
0 ⎟
.
⎟
⎠
0
1 n×n


We use the Kn −transformation that from Yn to Zn (= (zk )) by
Zn = K∗n Yn − Ȳ0 , K∗n = Pn C−1
n ,

(3.4)


where Ȳ0 = 1n · y0 ,

(3.5)

⎞
1
0 ··· 00
⎜ −1 10 ··· 0 ⎟⎟⎜−1⎜
0 ··· ⎟
Cn = ⎜ 0 −1 1
,
⎟
⎝ 0
⎠
0 −1 10000 −1 1
n×n
⎛

8

and
(3.6)

(n)
(n)
Pn = (pjk ) , pjk =


1
1
2
2π
(k − )(j − ) .
cos
2n + 1
2
2
n + 12
−1
Using the spectral decomposition C−1= Pn Dn Pn and Dn is a diagonal matrix
n Cn
2k−1
with the k-th element dk = 2[1 − cos(π( 2n+1 ))] (k = 1, · · · , n) , we write2k − 1∗2 π(k = 1, · · · , n) .
(3.7)
akn (= dk ) = 4 sin
2 2n + 1

Then the separating information maximum likelihood (SIML) estimator of Σx in(3.2) can be deﬁned byn
1 Gm = Σ̂x,SIM L =zk zk ,
mn k=1m(3.8)

where we set m = mn = [nα ] (0 < α < 1).
We need to use m terms in n and the reason for this becomes clear from the frequency domain analysis, which will be explained in Section 5.
Kunitomo, Awaya and Kurisu (2019) discussed the estimation of the variancecovariance matrix Σv when vi are i.i.d. vectors and some consistent estimatorsof Σv were developed. As we shall see in Section 5, the SIML estimation method(x)is quite robust, even when vi and vi are non-Gaussian stationary processes, andthey are serially-correlated.

3.2

General Non-Stationary Errors-in-Variables Model

In order to investigate non-stationary trend-cycle component, and stationary seasonality and measurement error component, we consider the general non-stationary
multivariate errors-in-variables model 4(3.9)

yi = xi + si + vi (i = 1, · · · , n) .

We take a positive integer s (s > 1), N , and n = sN for the resulting simplicity of
exposition and arguments, We explain the general model in three steps.
(i) The trend-cycle factor xi (i = 1, · · · , n) is a sequence of non-stationary I(1)process that satisﬁes(3.10)(x)Δxi = (1 − L)xi = vi

4,

It is possible to use the log-transformed data for multiplicative models. It is known thatthe standard model in X-12-ARIMA is a multiplicative one although it uses moving averages, forinstance.

9

with the lag-operator Lxi = xi−1 , Δ = 1 − L,
(x)
vi =

(3.11)∞
(x) (x)Cj ei−j ,j=0(x)(x) (x)

(x)

and ei is a sequence of i.i.d. random vectors with E(ei ) = 0 and E(ei ei ) =
(x)
(x)
Σ(x)
e (positive-semi-deﬁnite). The p×p coeﬃcient matrices Cj (= ckl (j)) are abso(x)(x)(x)
lutely summable and Cj  = O(ρj ), where 0 ≤ ρ < 1 and Cj  = maxk,l=1,··· ,p |ckl (j)|.
The random vectors vi (i = 1, · · · , n) are a sequence of stationary I(0) process withvi =(3.12)∞
(v) (v)

Cj ei−j ,

j=0
(v)

(v)

where the p × p coeﬃcient matrices Cj are absolutely summable and Cj  =(v)
O(ρj ), where 0 ≤ ρ < 1 and ei is a sequence of i.i.d. random vectors with(v)(v) (v)E(ei ) = 0, E(ei ei ) = Σ(v)e (positive deﬁnite).
(ii) The seasonal factor si (i = 1, · · · , n) is a sequence of stationary process 5 , which
satisﬁes(3.13)si =∞(s) (s)Csj ei−sj ,

j=0
(s)

where the lag-operator is deﬁned by Ls si = si−s (s ≥ 2), and ei is a a sequence(s)(s) (s)
of i.i.d. random vectors with E(ei ) = 0 and E(ei ei ) = Σ(s)e (a non-negative(s)deﬁnite matrix). The p × p coeﬃcient matrices Cj are absolutely summable and(s)Cj  = O(ρj ), where 0 ≤ ρ < 1.
(iii) Let fΔx (λ), fv (λ) and fs (λ) be the spectral density 6 (p × p) matrices of Δxi , vi
and si (i = 1, · · · , n) deﬁned by(3.14)fΔx (λ) = (

∞


(x)
Cj e2πiλj )Σ(x)
e (

j=0

(3.15)

fv (λ) = (

∞

j=0

(v)

Cj e2πiλj )Σ(v)
e (

∞


1
1
(x)
Cj e−2πiλj ) , (− ≤ λ ≤ ) ,
2
2
j=0

∞


1
1
(v)
Cj e−2πiλj ) , (− ≤ λ ≤ ) ,
2
2
j=0

5

When we assume that Δsi is stationary, most arguments in the following sections would go
through. We have omitted its details.
6
Our notation of spectral density is slightly diﬀerent from the standard notation used in Anderson (1971). Let μ = 2πλ and f A (μ) (−π ≤ μ ≤ π) be the spectral density in Chapter 7 ofAnderson (1971). Then, f (λ) = 2πf A (μ) (− 12 ≤ λ ≤ 12 ). The present deﬁtion of spectral density
corresponds to 2πf A (−μ) in some literature.

10

and
(3.16)

fs (λ) = (

∞


(s)
Csj e2πiλsj )Σ(s)
e (

j=0
(x)

(v)

∞


1
1
(s)
Csj e−2πiλsj ) (− ≤ λ ≤ ) ,
2
2
j=0

(s)

where we set C0 = C0 = C0 = Ip for normalization and i2 = −1 (see Chapter 7of Anderson (1971)).
Then the p×p spectral density matrix of the transformed vector process of diﬀerence
series Δyi (= yi − yi−1 ) can be represented as
fΔy (λ) = fΔx (λ) + (1 − e2πiλ )[fs (λ) + fv (λ)](1 − e−2πiλ ) .

(3.17)We denote the long-run variance-covariance matrices of trend-cycle and noise components for g, h = 1, · · · , p as
(x)

(v)

Σx = fΔx (0) (= (σgh )) , Σv = fv (0) = (σgh ) .

(3.18)
respectively.

One important often neglected issue is that when applying the diﬀerencing procedure to non-stationary time series and using the standard statistical method for
multivariate stationary time series, there is no guarantee of keeping the relationships
among the original time series as they were by using the transformations. Although
Engle and Granger (1987) and Johansen (1995) noticed this problem, they did notconsider the frequency domain aspect with seasonality and measurement errors. The
SIML ﬁltering approach may shed a new light on the relationships among the time
domain and frequency decompositions of non-stationary multivariate time series.

4

The SIML Filtering Method

4.1

Basic Filtering

We introduce the general ﬁltering procedure based on the Kn −transformation in(3.4). When we interpret that the elements of the resulting n × p random matrix
Zn take real values in the frequency domain, it is easy to understand their roles.
Since Pn is a kind of real-valued discrete Fourier transformation, vectors zk (k =1, · · · , n) in Zn are asymptotically uncorrelated, as we shall discuss in Section 5.1.
We investigate the partial inversion of the transformed orthogonal processes. Let
an n × p matrix(4.19)X̂n (Q) = Cn Pn Qn Pn C−1
n (Yn − Ȳ0 )

11

and
(4.20)∗
Zn = Pn C−1
n (Yn − Ȳ0 ) , Yn = Ȳ0 + Xn + Sn + Vn ,where X∗n = (x∗t ), Sn = (st ) and Vn = (vt ) are n × p matrices, x∗t = xt − x0 (t =1, · · · , n) and x0 = y0 is the initial vector. (We use the subscript t instead of i inthis subsection.)
The stochastic process Zn is the orthogonal decomposition of the original time series
Yn in the frequency domain and Qn is an n × n ﬁltering matrix. Because Yn consistof non-stationary time series, we need a special form of transformation Kn in (3.4).
We give explicit forms of two examples, including the trend-cycle ﬁltering and theband ﬁltering procedures. Although there can be many possible ﬁltering procedures
within our general framework, it is useful to discuss some linear ﬁltering procedures.
Let an n × n diagonal matrix
Qn =(4.21)n(n) (n)

wt,n et et

t=1
(n)



and et = (0, · · · , 1, · · · ) (t = 1, · · · , n) are the unit vectors and wt,n (t = 1, · · · , n)are some non-negative constants.
We start with the case when wt,n = 1 (t = 1, · · · , n) and we have the identity
matrix Qn = In . Then we ﬁnd Cn Pn Qn Pn C−1
n = In . There can be useful cases andwe present two cases as the trend-cycle ﬁltering and the band ﬁltering by choosing
wt,n = 1 or 0 for some t s. Although the ﬁrst example could be regarded a special
case of the second one, the analysis of trend-cycle component in the ﬁrst example
has an important role in economic time series.
(i) Trend-Cycle Filtering: Let an m × n (m < n) choice matrix Jm = (Im , O),
and let also n × p matrix(4.22)

= Cn Pn Jm Jm Pn C−1X̂(m)nn (Yn − Ȳ0 )

and an n × n matrix Qn = Jm Jm .
We construct an estimator of n×p hidden state matrix X∗n only in the lower frequency
parts by using the inverse transformation of Zn and deleting the estimated seasonaland noise parts. We denote the hidden trend-cycle state based on m frequencies as(4.23)

∗= Cn Pn Jm Jm Pn C−1
X(m)nn Xn .

This quantity is diﬀerent from X∗n because xt (t = 1, · · · , n) in (3.9) and (3.10)
contains not only the trend-cycle component of yt (t = 1, · · · , n), but also the noise
component in the frequency domain, which is diﬀerent from the measurement noise
component vt (t = 1, · · · , n) in (3.9)-(3.13). We try to estimate the trend-cycle
component of xt by using (4.22) and recover the trend-cycle component of Xn near12

at the zero frequency because the eﬀects of diﬀerenced measurement error noises(vt − vt−1 ) are negligible around at zero frequency. This method diﬀers from someexisting procedures that consider the decomposition of time series only in the time
domain. Our arguments can be justiﬁed by using the frequency decomposition of yt(n)and rt = Δyt (= yt − yt−1 and y0 being ﬁxed), and we shall discuss this issue inSection 5.2.
We partition Pn into [m + (n − m)] × [m + (n − m)] matrices as
P11 P12Pn =P21 P22and thenP n Jm J m P n =

P11
P21




(P11 , P12 ) = In −

P12
P22


(P21 , P22 ) .
After straightforward calculations (see the Appendix for the derivation), the (j, j )(n,m)
th element of An = Pn Jm Jm Pn (= (aj,j  )) is given by(n,m)(4.24) aj,j(n,m)

aj,j 


2mπ(2j − 1)
sin 2n+1
2m
1
,
=
+
π
2n + 1 2n + 1 sin 2n+1
(2j − 1)




2mπ
2mπ
(j + j − 1) sin 2n+1
(j − j )sin 2n+11+(j = j ) .
=


π
π
2n + 1 sin 2n+1 (j + j − 1) sin 2n+1 (j − j )(m)

It is possible to evaluate MSE of statistical state vector estimation. Let X̂n =
(m)
(m)
(m)
(m)
(m)
(X̂ni ) (X̂ni is an n × 1 vector, i = 1, · · · , p) and Xn = (Xni ) (Xni is ann × 1 vector for i = 1, · · · , p). By decomposing Yn − Ȳ0 = X∗n + Sn + Vn and(m)X̂n − X∗n = Cn Pn Jm Jm Pn C−1
n (Sn + Vn ), we ﬁnd
(m)

(m)

(m)

(m) 





∗ (s+v)
E[(X̂ni − Xni )(X̂nj − Xnj ) ] = K∗−1
(i, j)K∗n Qn K∗n −1 ,n Qn Kn Γn(s+v)
where K∗n = Pn C−1(i, j) is the n × n variance-covariance matrix of Sni + Vni
n , Γnand Snj + Vnj (n × 1 vectors for i, j = 1, · · · , p), and, Sni and Vni (n × 1 vectors)
are the i−th and j−the column vectors of Sn and Vn , respectively.
Since (25) does not depend on xt , (22) minimizes the MSE with respect to unknown
state vector xt to estimate (23) and it is optimal in this sense.

(ii) Band Filtering : We consider a general ﬁltering based on the Kn transformation in (3.4) and use the inversion of some frequency parts of the random matrix
Zn . The leading example is the analysis of seasonal frequencies in the discrete timeseries and we take s (> 1) being a positive integer as the seasonal lag.
13

Let an m2 × [m1 + m2 + (n − m1 − m2 )] choice matrix Jm1 ,m2 = (O, Im2 , O) (we takem1 + m2 < n), and let also n × p matrix


1 ,m2 )= Cn Pn Jm1 ,m2 Jm1 ,m2 Pn C−1X̂(mnn (Yn − Ȳ0 )(4.25)

and an n × n matrix Qn = Jm1 ,m2 Jm1 ,m2 .
When we have a particular seasonal frequency s (> 1), for instance, we can takem1 = [2n/s] − [m/2] and m2 = m. We set s = 4 for quarterly data and s = 12 formonthly data. (We may have many seasonal frequencies in data analysis as pointed
out by Granger and Hatanaka (1964) already.)
As in the trend-cycle ﬁltering problem, (4.25) is the SIML-ﬁltering value for∗
1 ,m2 )X(m= Cn Pn Jm1 ,m2 Jm1 ,m2 Pn C−1nn (Xn + Sn ) ,(4.26)and it is an estimate of some frequency components of xt + st (t = 1, · · · , n) in (3.9),(3.10) and (3.13). The ﬁltering problem becomes diﬃcult because some frequency
component of yt includes not only some component xt at the same frequency, butalso some component of the measurement errors at the same frequency.
After straightforward calculations (see the Appendix) in this case, the (j, j )-th(n,m ,m )element of An = Pn Jm1 ,m2 Jm1 ,m2 Pn (= (aj,j  1 2 )) is given by(n,m1 ,m2 )(4.27)aj,j(n,m ,m2 )aj,j  1
1 +m2 )π1 )π(2j − 1) − sin 2(m(2j−1)
sin 2(m2n+1
1
2m2
2n+1
+
,
=
π
2n + 1 2n + 1
sin 2n+1(2j − 1)

1 +m2 )π1 )π(j + j − 1) − sin 2(m(j + j − 1)
sin 2(m2n+1
1
2n+1
=
π
2n + 1
sin 2n+1(j + j  − 1)

1 +m2 )π1 )πsin 2(m2n+1(j − j ) − sin 2(m
(j
−
j
)

2n+1
+
(j = j ) .
π
sin 2n+1 (j − j )When m1 = 0 and m2 = m, the resulting formula reduces to the trend-cycle ﬁlteringcase. There are existing ﬁltering procedures having the frequency domain interpretation, but it seems that our procedure diﬀers from some existing literature becauseof (4.24) and (4.27).
It is also possible to evaluate MSE of the state vector estimation in the same way(m ,m )as Case (i). Let X̂n 1 2 be an estimate of the hidden state as(s,m ,m )∗Xn 1 2 = Cn Pn Jm1 ,m2 ,n Jm1 ,m2 ,n Pn C−1
n (Xn + Sn ).
By using the similar calculation as Case (i) and the notation of Γ(v)n (i, j) (i, j =1, · · · , p), we ﬁnd(4.28)(m ,m2 )E[(X̂ni 1(s,m1 ,m2 )

− Xni



(m ,m2 ))(X̂nj 1


∗ (v)
∗
∗ −1
,
= K∗−1n Qn Kn Γn (i, j)Kn Qn Kn

14(s,m1 ,m2 )

− Xnj



)) ]

where Γ(v)n (i, j) is the n × n variance-covariance matrix of Vni and Vnj . We use an(m2 − m1 ) × n choice matrix Jm1 ,m2 = (O, Im2 −m1 , O) and Q∗n = In − Jm1 ,m2 Jm1 ,m2 .
If there are several seasonal frequencies, we need more complicated ﬁltering procedures.

5

On a Statistical Foundation

5.1

Underlying Asymptotic Theory

At ﬁrst glance, the SIML ﬁltering procedure may be regarded as an ad-hoc statisticalprocedure without any mathematical foundation. However, it has a rather solid
statistical foundation.

(n)
2π
1
(j − 12 )(k − 12 ) , pjk = √2n+1(eiθjk + e−iθjk ) and for Yn = (yi ) weLet θjk = 2n+1
write zk (k = 1, · · · , n) as(5.29)(n)zn (λk ) =n(n)
pjk rj , rj = yj − yj−1 ,j=1which is a (real-valued) Fourier-transformation and y0 is ﬁxed.
(n)Then, we ﬁnd that zn (λk ) (k = 1, · · · , n) are the (real-valued) Fourier-transformation(n)of data at the frequency λk (= (k − 1/2)/(2n + 1)), which is a (real-part of) estimate of the orthogonal incremental process z(λ) (0 ≤ λ ≤ 1/2), which is continuous
in the frequency domain.
We shall utilize an asymptotic theory for the stationary linear processes because
the time series model deﬁned by (3.9)-(3.13) can be regarded as a special case. Let(5.30)

ri = Δyi = μ +∞Γj wi−j ,j=0where μ is a constant vector (we assume μ = 0 for simplicity in this section), Γjare p × p matrices, and wi are a sequence of mutually independent random variables
with E[wi ] = 0, E[wi wi ] = Σu (> 0). The errors-in-variables model in Section 3
implies that the p × p matrices Γj satisfy ∞h=0 Γj  < ∞ .
(n)We then summarize the useful result on z(λk ). Although it could be regarded asa direct extension of Theorem 8.4.3 of Anderson (1971) for discrete and (ergodic)
stationary time series, we could not ﬁnd the following representation. The proof isgiven in the Appendix.
Proposition 1: Letrj (j = 1, · · · , n) be an ergodic stationary stochastic process
given by (5.30) with ∞h=0 Γj  < ∞ and the fourth order moments of each elementof wi are ﬁnite.
15

n
(n)
(n)
Let also zn (λk ) =j=1 pjk rj and rj be an ergodic stationary sequence withE[rj ] = 0, Γ(h) = E(rj rj−h ), and the (symmetrized real-valued) spectral density
matrix(5.31)fSR (λ) =∞cos(2πhλ)Γ(h) , 0 ≤ λ ≤

h=−∞

1
2

is the positive deﬁnite and bounded (real-valued and symmetrized) spectral matrix.
(n)(n)
Assume that λk → s, λk → t as n → ∞ for 0 < s < t < 12 . Then, as n −→ ∞



(n)
zn (λk )0
fSR (s)w(5.32)−→ N2p 0, [
] .
(n)
0
fSR (t)zn (λk )This proposition covers the general model with (3.9)-(3.13) with the moment conditions because Δyi are stationary. As the asymptotic variance-covariance matrix(n)of the orthogonal random vectors zn (λk ) is the (symmetrized real) spectral density
matrix, it can be estimated consistently. When we have noise terms as in Section 3,
it is not possible to estimate the (long-run) variance-covariance matrix Σx of trendcycle component simply by using the diﬀerenced time series rj (= Δyj ) = yj − yj−1 .
It is because(5.33)



E[rj rj ] = E[Δxj Δxj ] + E[Δsj Δsj ] + E[Δvj Δvj ] .

The SIML estimator of fSR (0) (= fΔx (0)) = Σx in the general case can be deﬁned
by
1 
(n)
(n) 
Gm =
(zn (λk ))(zn (λk )) .
mn k=1m(5.34)

Then we summarize the basic property of the SIML estimation of Σx and the derivation is given in the Appendix.
(x)

Proposition 2: Assume that the fourth order moments of each element of vi , vi
(s)and vi (i = 1, · · · , n) in (3.9)-(3.13) are bounded. We set mn = [nα ] (0 < α < 1)).
Then, as n −→ ∞(5.35)

p

Gm −→ Σx .

The above result has some implication on the use of the frequency domain analysisof (non-stationary) multiple time series. For 0 ≤ μ ≤ 12 , let the symmetrized (realvalued) spectral density matrix of Δyj (j = 1, · · · , n) be fSR,Δy (μ) = (1/2)[fΔy (μ)+
16f¯Δy (μ)] , where f¯( · ) is the complex conjugate of f and we have the initial conditiony0 .
Similarly, we ﬁnd that fSR,Δx (λ) = (1/2)[fΔx (λ) + f¯Δx (λ)] ,fSR,Δs (λ) = (1/2)(1 − e2πiλ )[fs (λ) + f¯s (λ)](1 − e−2πiλ ) andfSR,Δv (λ) = (1/2)(1 − e2πiλ )[fv (λ) + f¯v (λ)](1 − e−2πiλ ).
They are the consequence of (3.17), (5.29) and (5.31) for the symmetrized spectral(real-valued) density matrices.
From this interpretation, it may be easy to ﬁnd that Gm is a consistent estimatorof the long-run variance-covariance matrix, which is the spectral density matrix at
the zero-frequency fSR,Δx (0) when we use the Zn − transformation of data.

5.2

Frequency Interpretation of SIML-ﬁltering

In the traditional statistical time series analysis for a stationary discrete (vector)
process with the (complex-valued) spectral distribution F , it has a representation
with right-continuous (complex-valued) orthogonal increments in the frequency domain (see Doob (1953), Brockwell and Davis (1990) for the details). Chapter 7 ofAnderson (1971) is informative because of its discussion on real-valued representations although it has only univariate cases. The real-valued multivariate orthogonal
processes and the spectral density matrix play important roles in our formulation.
(n)For λk = (k − 1/2)/(2n + 1) (k = 1, · · · , n), we rewrite (5.29) as(5.36)(n)zn (λk ) =

n

1
2
(n)
(n)
rj √
cos[2πλk (j − )] (k = 1, · · · , n).
22n + 1j=1

Then, by using the inversion transformation with Pn , we can conﬁrm that
(5.37)

r(n)
s =

n


(n)

psk zn (λk ) (s = 1, · · · , n).

k=1
(n)

It is another representation of Rn = (ri ) = C−1n X̂n (Q) in (4.19) when Qn = In .
(n)
For any s (s = 1, · · · , n), rs can be recovered as the weighted sum of othogonal(n)(n)
processes zn (λk ) at frequency λk (k = 1, · · · , n). We then, by using Yn = Cn Rn ,(n)recover the non-stationary process yt (t = 1, · · · , n) given the initial condition y0
as(5.38)(n)yt= y0 +

t


r(n)
.
s

s=1

Let
(5.39)11
αn (λ(n)m ,j − ) =
2

m


1
(n)
[2 cos 2πλk (j − )]
n k=1
2
17

(n)

Then, when λm → λ as n → ∞ (0 < λ < 12 ), using Lemma 5.1 of Kunitomo et al.
(2018), we ﬁnd
2 sin 2πλ(j − 12 )
1
1
,
j
−
.
αn (λ(n)
)
→
α(λ,
j
−
)
=
m
2
2
π(j − 12 )If we set the uncorrelated stochastic process of uncorrelated increments with con(n)1tinuous parameter λ (0 ≤ λ ≤ 2 ) as An (λ) = nj=1 α(λ, j − 12 )rj , then we ﬁnd 1
2(5.40)
01
cos[2πλ(s − )]dAn (λ) = r(n)(s = 1, · · · , n) .
s
2This corresponds to the continuous representation of a discrete (real-valued) stationary time series in the frequency domain (see Chapter 7.4 of Anderson (1971)).
If we write the limit of A(λ) = limn→∞ An (λ) (assuming it exists), the (real-valued)
spectral distribution matrix FRS for any 0 ≤ λ1 < λ2 ≤ 1/2 can be deﬁned as λ2

(5.41)FRS (λ2 − λ1 ) = E[(A(λ2 − λ1 )A(λ2 − λ1 ) ] =fRS (λ)dλλ1

if FRS is absolutely continuous and the matrix-valued density process fRS (λ) (0 ≤λ1 < λ2 ≤ 1/2) exists.
(m,n)
(m,n)) = C−1
are p × 1 vectors for
From (4.22), we set R̂n (m) = (r̂in X̂n (m) and r̂ii = 1, · · · , n. If we write=r̂(m,n)s(5.42)m(n)

psk zn (λk ) (s = 1, · · · , m; 0 < m < n),

k=1
(m,n)

it is the trend-cycle estimate of the SIML-ﬁltering value for rs
corresponding element of Rn (m) (= Cn Xn (m)). It corresponds tor(m,n)=
s(5.43)m, which is the(n)

psk z∗n (λk ) (s = 1, · · · , m; 0 < m < n),

k=1
(n)

where z∗n (λk ) are constructed from the n × p hidden states matrix Xn instead of
the observed n × p matrix data Yn . Hence (5.42) is the same as the element of
(m)
(n)
C−1
in (4.22), and for λm = m/(2n) in the frequency domain it is a discreten X̂nversion of λ(n)m1(n) (n)(5.44)r̂s (λm ) =cos[2πλ(s − )]dAn (λ) .
20(m)in (4.23), and it has the correThen, (5.43) is the same as the element of C−1
n Xn
sponding (continuous) version in the frequency domain.
18

m1 +m2(m ,m ,n)(n)(m2 ,n)(m ,n)
Similarly, r̂s 1 2=− r̂s 1
k=m1 +1 psk zn (λk ) = r̂s
0 < m1 < m2 < n) can be regarded as a discrete version of λ(n)m2
1(n) (n)(n)cos[2πλ(s − )]dAn (λ) .
(5.45)r̂s (λm1 , λm2 ) =
(n)
2
λ m1

(s = 1, · · · , m;

From our interpretation of the SIML ﬁltering we ﬁnd an interesting representation of (discrete time and real-valued) stationary processes and orthogonal incremental stochastic processes. They are closely related to our method of data analysis
for non-stationary vector time series.

6

Some Simulation and the Choice of Frequencies

6.1

A Guide of Choosing Frequencies

When we are interested in ﬁltering a non-stationary time series with trend-cycle,
seasonality and measurement errors, we need to choose the parameter m. We give a
guide to set m in the trend-cycle ﬁltering case for practical purpose. From the dis(n)cussion of Section 5.2, the orthogonal process zn (λk ) corresponds to the frequency(n)λk = (k − 12 )/(2n + 1) (k = 1, · · · , n). When we are interested in the trend-cycle(n)component, we may only use λk = (k − 12 )/(2n + 1) (k = 1, · · · , m) and then the
maximum frequency is approximatelym.
(6.46)
λ(n)
max =
2n
For instance, when we have monthly data over 20 years as an example, we haven = 240 and s = 12. Since we have seasonal frequencies, we want to ﬁnd trendcycle components as business cycles more than 1.5 year, say. Then an appropriate(n)maximum frequency would be λmax = 1.5/24 and then we could take m∗ = 480 ×(1.5/24) = 30.
As another example if we have quarterly data over 30 years, we have n = 120 ands = 4. Since we have seasonal frequency, we want to ﬁnd trend-cycle components as
business cycles more than 1.5 year, say. Then an appropriate maximum frequency(n)would be λmax = 1.5/8 and then we could take m∗ = 240 × (1.5/8) = 45.
If we were interested in the trend-cycle component of the non-stationary time series,
these choices might be reasonable candidates.
As we shall see in the next subsection, this point could be checked by simulation
for prediction.

6.2

Some Simulation(m)

When we have estimates of the state variables xi (i = 1, · · · , n), the estimates(m)(m)
of error components are v̂i = yi − x̂i (i = 1, . . . , n). Then, an estimated MSE
19

of the one-step ahead prediction errors based on the SIML-smoothing or ﬁltering isgiven by

PMSEn (h) = E[(yn+h − yn )(yn+h − yn ) |Fn ] ,(6.47)

where Fn is the σ−ﬁeld (information) available at n.
One may try to minimize the estimated h-step prediction MSE by choosing anappropriate m. It may be reasonable to take h = 2, 3 as short-run prediction andh = 4, 8 for long-run prediction from our limited experiments.
We only show a result on the simple trend plus noise model when p = 1, xi =(x)(x)
xi−1 + vi and yi = xi + vi (i = 1, · · · , n) with vi ∼ N (0, σx2 ) and vi ∼ N (0, σv2 ).
The criterion function is the prediction MSE given by
1
PMSE∗n (h) =

n


(m)

h i=n−h+1(m)(yi − x̂i )2 ,(m)where x̂i is the estimate of xi .
We present the minimum m as m∗ based on the trend-cycle ﬁltering as Tables 6.1-6.3
by taking h = 2, 4, 6, 8 and n = 80, 120, 200, 400. In our simulations, as n increases,
we have larger choice of m∗ while m∗ decreases as h increases. When we have a longhorizon with h, it may be natural to use a small number of lower frequencies, andfor a wide range of σx and σv .
Table 6.1 : PMSE and Choice of m
(σx = 0.3 , σv = 0.05)
n
h=2
h=3
h=4
h=5
h=6
h=7
h=8

80
12
8
6
5
4
4
3

120
19
12
9
7
6
5
4

20

200
32
21
15
13
10
9
8

400
65
42
32
26
21
18
16

Table 6.2 : PMSE and Choice of m
(σx = 0.3 , σv = 0.4)
n
h=2
h=3
h=4
h=5
h=6
h=7
h=8

80
13
8
6
5
4
4
3

120
19
13
9
8
6
5
5

200
33
21
16
13
11
9
8

400
66
43
32
26
21
18
16

Table 6.3 : PMSE and Choice of m
(σx = 0.3 , σv = 1.0)
n
h=2
h=3
h=4
h=5
h=6
h=7
h=8

6.3

80
13
9
6
5
4
4
4

120
20
13
10
8
7
6
5

200
33
22
16
13
11
9
8

400
67
44
33
26
22
19
16

A Comparison of the SIML Filtering and HP Filtering

In order to characterize some property of the SIML ﬁltering, we give an illustrative
example. For this purpose, we take a monthly original consumption data in Section 2(Kakei-Chosa series in Figure 3) and show the result of the SIML ﬁltering and the HP
ﬁltering with two parameter values (λ) as Figure 4. (We used hpf ilter procedurein mFilter R-package for HP ﬁlter and took m = 12 for the SIML ﬁlter, which(n)corresponds to frequencies over 2 year cycles since λmax = 1/24.) An important pointis that when we have strong seasonality in data, the SIML often gives reasonableestimates while the usefulness of HP ﬁlter depends on the the particular parametervalue chosen for λ.
Further comparison would be interesting, but it is beyond the scope of this paper.

7

Applications

7.1

An interpretation of Müller and Watson (2018)

Recently, Müller and Watson (2018) have proposed the so-called long-run co-variabilityof macro-economic time series. They have investigated many non-stationary time21

Figure 4: An Example of SIML and HP Filters. (An illustrated comparison of the
SIML and HP Filters by using consumption data.)

22

series and obtained some empirical ﬁndings. In this sub-section, an interpretationof their method will be given as measuring the relationships among long-run trendcycle in our framework when p = 2. Moreover, we obtain an important consequence
from Proposition 2 in Section 5.1.
In this sub-section, we consider the case of p = 2 in (3.9)-(3.13). Let 2×2 matrices
(x)
(x)
(x)
Σx = (σij ), we then deﬁne the (long-run) regression coeﬃcient β = [σ22 ]−1 σ21(x)(m)
under the assumption that σ22 > 0. Let also Gm (0) = (gij ), and an n × 2 matrix(w1 , w2 ) = C−1
n (Yn − Y0 ) .

(7.48)For estimating β, we deﬁne(7.49)

(m)

(m)









β̂ = [g22 ]−1 g21 = [w2 Pn Jm Jm Pn w2 ]−1 [w2 Pn Jm Jm Pn w1 ] .

This quantity can be interpreted as the least squares slope of the transformed vector
from y1n on the transformed vector from y2n for an n × 2 matrix Yn = (y1n , y2n ),which is essentially the same as the one proposed by Müller and Watson (2018) 7 .
Then, from Proposition 2 in Section 5 and its proof in the Appendix we immediately
obtain the following result.
Proposition 3 : In (3.9)-(3.13) with p = 2, we assume that the fourth order(x)(s)
moments of each element of vi , vi and vi (i = 1, · · · , n) are bounded, and Σx is
positive deﬁnite.
(i) We ﬁx an m, then β̂ is not consistent when n → ∞.
(ii) Set mn = [nα ] and 0 < α < 1, then as n −→ ∞pβ̂ − β −→ 0 .

(7.50)The second part of this proposition is an extended version of the ﬁrst part of Theorem
4.1 of Kunitomo and Sato (2017).
A Simulation Example
To illustrate our arguments in Proposition 3, we performed a set of Monte Carloexperiments under the simple situation with p = 2 (the replication was 10,00 in eachsimulation). The model is given by(x)(7.51)

xi = (vx1,ix) = ( 1,i−1 ) + ( 1,i(x) ) ,x2,ix2,i−1v2,i

7

In their notation, m corresponds to q, which is ﬁxed. They did use the (diﬀerenced) stationary
data, and thus, we could interpret that they calculated the linear regression from the ﬁltered dataX̂∗n = Pn Jm Jm Pn C−1
n (Yn − Y0 ) as a modiﬁcation of (4.19) in our notation.

23

(7.52)yi = (y1,ix + βx2,iv) = ( 1,i
) + ( 1,i ) ,
y2,i
x2,i
v2,i
(x)

(x)

(x)



where we have generated the normal errors vi = (v1,i , v2,i ), vi = (v1,i , v2,i ), (i =(x)(v)
1, · · · , n) with zero means, variances σj = 0.75 and σj = 1.5 (j = 1, 2), zerocovariances and an initial value x0 .
This is a two dimensional I(1)-process, which is not co-integrated, but close to a
cointegrated-process. We present the ﬁnite sample properties of the (naive) leastsquares (LS) estimator from the original data, the Müller-Watson (MW) estimatot,
and the SIML (SIML) estimator from the transformed data.
The simulation results in Tables 7.1-7.3 are consistent with our arguments in
Proposition 3. In the tables, Direct LS stands for the standard regression on y1n on
y2n , while MW and SIML stand for Muller-Watson and SIML, respectively. When
we have two-dimensional time series and they are not cointegrated, the standard
least square method is badly biased. The procedure proposed by Muller and Watson(2018) often gives reasonable results, but its variance does not decrease as the sample
size increases. The SIML estimation method has reasonable ﬁnite sample properties
as well as reasonable asymptotic properties and it is applicable to more generalcases.
Table 7.1 : Simulations (Direct LS)
n
true
mean
SD

200
1000
1.2000 1.2000
2.2016 2.2287
0.0930 0.2183

30000
1.20002.09181.4901

Table 7.2 : Simulations (MW)
n
true
mean
SD

200
1000
1.2000 1.2000
1.1418 1.2129
0.5279 0.5301

30000
1.20001.18750.5337

Table 7.3 : Simulations (SIML)
n
true
mean
SD

200
1000
1.2000 1.2000
1.0828 1.1807
0.3278 0.1904

24

30000
1.20001.20200.0766

7.2

Some Macro-economic data and Constructing a Monthly
Consumption Index

As illustrations for the empirical analysis, we have used our ﬁltering method toanalyze Japanese quarterly (real) consumption-GDP data as the ﬁrst example, andthree sets of monthly consumption data as the second example, which have been
discussed in Section 2. Appendix B contains all ﬁgures in this subsection.
As the ﬁrst example, the ratio of real GDP and real consumption (quarterly
original series) and their time series plots are given in Figures 1 and 2. They show
non-stationarity in their trend-cycle and seasonality as typical macro-economic variables. We then calculate the transformed data using Pn (a kind of real Fouriertransformation) as Figure 5 from the original consumption data. In this case, thetransformed series gives a wild up and down ﬂuctuation, and trend-cycle and seasonality are sucked. To ﬁnd seasonal components from the original series, we calculated
the realized zk (k = 1, · · · , n) (Figure 6) and the empirical cumulative distribution of
zk2 (k = 1, · · · , n) (Figure 7), which roughtly correspond to the normalized and realvalued sample spectral distribution function. Each component of zk (k = 1, · · · , n)
is an orthogonal decomposition in the frequency domain. Because we have quarterly(n)macro-data, we have a large up and down around (λk =) 0.25, which corresponds tothe seasonal frequency at s = 4. The empirical spectral distribution has an abruptchange at this frequency. From these ﬁgures, we can judge that the real Fouriertransformation based on Kn does give useful information.
As it has been a practice in time series data analysis to use seasonal diﬀerencing
Δs yi = yi − yi−s (s = 4) in the Box-Jenkins method, we calculated the real Fouriertransformation based on Pn (Figure 8) after seasonal diﬀerencing. Although thecontribution of the resulting orthogonal process around the seasonal frequency could
be signiﬁcant, there are some rather wild ﬂuctuations on at many other frequencies
by using Pn −transformation. Because we have some diﬃculty in interpreting theresulting time series, it may not be possible to justify the seasonal diﬀerencing
procedure, and we recommend not to use this representation. In the following
analysis, we simply use the diﬀerencing and then use the frequency domain analysis.
In Figure 9, we have analyzed real Quarterly-GDP. We show one example withm = [n.99 ] and the deleted seasonal frequency are around 48-52 (48/196-52/196 in[0, 1/2]) and we delete extremely high-frequency part. We deleted only ﬁve transformed data around the seasonal frequency and the main intention was to investigate
the eﬀect of seasonality with a narrow band. We have taken α = .99 since we wanted
to remove some contribution of high frequency, but we could have used other choices
and the results are not much diﬀerent from Figure 8 in most cases. We compared theﬁltered time series using our method and the oﬃcial (published) seasonally adjusted
time series. We found that the diﬀerences in two time series (i.e., the publishedtime series and the SIML ﬁltered time series) are rather small and they are often
of negligible magnitude. Although our ﬁltering procedure is simple, this empirical25

example suggests the usefulness of our method developed in this study.
As the second example, we have analyzed three consumption (monthly) timeseries and the quarterly consumption time series, which are mentioned in Section 2.
Then this is an empirical example when p = 4 and s = 12 with missing observations
because the quarterly one series cannot give monthly information. As we have seenin Figure 3, three monthly consumption series have similarities and some diﬀerencesin their components. In our example, our goal is to construct the monthly consumption index, which is consistent with the observed quarterly consumption timeseries in trend-cycle component. Because of non-stationary trend-cycle, seasonaland measurement errors, it may not be obvious to construct a useful consumption
index using existing statistical tools.
Let Yi (i = 1, · · · , n) be the target (quarterly) time series and Zjt (j = 1, 2, 3; t =3(i − 1) + l, l = 1, 2, 3) be the j-th monthly time series (t = 0 is the initial period
and we ﬁx the initial values Y0 , Zj0 ). Then the criterion function is(7.53)

2
n
3
3



(T )
(T )
ΔŶi −wj (ΔZj,3(i−1)+l ) ,MSE(m, m1 , m2 , m3 , w1 , w2 , w3 ) =
i=1
(T )

(T )

(T )

j=1

l=1

where ΔŶi
= Ŷi − Ŷi−1 , (the trend part of the estimated ΔYi because we
(T )
(T )
(T )
observe quarterly data on Yi ), ΔZjt = Zjt − Zj,t−1 (the trend parts of ΔZjt ), andwj (j = 1, 2, 3) are (unknown) weight coeﬃcients and m, mj (j = 1, 2, 3) are thenumbers of trend-cycle ﬁltering. In the above formulation, we need to measure theprediction errors based on diﬀerenced data because we have non-stationary trendcycle component.
Using the least squares method, we minimized the MSE criterion with respect
to the underlying parameters. The estimated wj (j = 1, 2, 3) are 3.69, 5.19, and1.64 (while the measurement units are diﬀerent), but their magnitudes are comparable to the published quarterly consumption level at 2002Q1-2016Q4), which arestatistically signiﬁcant at 1%. We have chosen m = 29, m1 = 36, m2 = 23 andm3 = 33. Although it may be possible to use other possibilities, but in our limited
experiments, we found some improvements in prediction error over other cases with
diﬀerent combinations of m and mj (j = 1, 2, 3).
The black curves are the original series and the red curves are the estimatedtrend curves in Figures 10 and 11 for two monthly series. By taking relatively large
mj (j = 1, 2, 3), we can recover the cycle components of each series, which are
crucial as the indicators of macro-business condition. In Figure 12, the green curveshows the predicted state value calculated from the latest observed (quarterly) dataplus the predicted monthly part based on the estimated parameters. As there is
no monthly observation of quarterly published consumption, we draw their latest(quarterly) level with the black curve and the estimated SIML (ﬁltered) values with26

the red curve 8 . Overall, we found that while our procedure is much simpler than the
X-12-ARIMA seasonal adjustment with reg-ARIMA model, the estimated ﬁlteredseries are consistent with the (published) quarterly series given the fact we have nonstationary trend-cycle, seasonal and measurement errors components even when we
do not have large samples. In Figure 13. we have drawn the prediction errors interms of the diﬀerenced value Yi (i = 1, · · · , n) based on our procedure. This ﬁgure
illustrates the usefulness of the procedure because the macro-economic time series
are non-stationary with measurement errors.
The empirical data analyses in this subsection are presented for illustrations. We
are currently investigating the consumption data and the results will be reported onanother occasion.

8

Concluding Remarks

When the observed non-stationary multivariate time series contain seasonality andnoises, it may be diﬃcult to disentangle the eﬀects of trend-cycle, seasonal and measurement errors. In real (seasonally unadjusted) times series, we often observe nonstationary trend-cycle, seasonality and measurement errors while the X-12-ARIMA
program in oﬃcial agencies uses the univariate reg-ARIMA model to remove theseasonality from original time series. In this study we investigate a new procedure
to decompose time series into non-stationary trend-cycle components, stationary
seasonal and noise (or measurement errors) components. The resulting method fornon-stationary multivariate series is simple and free from the underlying distributions of components. Hence, it is robust against possible misspeciﬁcation in the
non-stationary multivariate economic time series. An important conclusion is that
it is useful to transform the observed time series using the Kn −transformation andinvestigate the transformed Zn series.
In empirical example in Section 7, we have illustrated our method to analyze
quarterly and monthly macro-consumption data in Japan. We presented a way toconstruct the monthly consumption index as the second example, which is consistent
with the published or oﬃcial (GDP-)consumption quarterly data. Although the
problem is practically complicated, we have shown that our method gives a usefulway for practical purpose.
There can be several further problems. Although it is easy to handle the Kn −transformations of non-stationary multivariate time series and construct the transformed Zn −data, there is a non-trivial initial value problem of ﬁltering. Anotherdirection would be to handle some abrupt changes in multiple time series, such as
the changes of consumption tax in Japan. Some progress on these problem has been
8

One notable event was the introduction of consumption tax in April 2014 and a sharp deviationof trend. In the present study, we do not focus on this aspect and it will be in Sato and Kunitomo(2020b).

27

made (Sato and Kunitomo (2020a,b)), which will be reported in other occasions. As
there are many important empirical applications, we need to develop a systematic
statistical procedure.

References[1] Anderson, T.W. (1971), The Statistical Analysis of Time Series, Wiley.
[2] Brillinger, D. and M. Hatanaka (1969) ”An Harmonic Analysis of Nonstationary
Multivariate Economic Processes,” Econometrica, Vol. 35, 131-141.
[3] Brillinger, D. (1980) Time Series : data analysis and theory, Expanded Edition,
Holden-Day.
[4] Brockwell, P. and R. Davis (1990), Time Series : Theory and Methods, 2nd
Edition, Wiley.
[5] Baxter, H. and R. King, (1999) ”Measuring Business Cycles : Approximate
Band-Pass Filters for Economic Time Series,” Review of Economics and Statistics, Vol.81-4, 575-593.
[6] Christiano, L. and T. Fitzgerald (2003) ”The Band Pass Filter,” International
Economic Review, Vol.44-2, 435-465.
[7] Doob, J.L (1953), Stochastic Processes, Wiley.
[8] Durrett, R. (1991), Probability : Theory and Examples, Duxbury Press.
[9] Engle, R. and C.W.J. Granger (1987), “Co-integration and Error Correction,”
Econometrica, Vol.55, 251-276.
[10] Findley, D., Monsell, B., Bell, W., Otto, M. and Chen, B. (1998),”New Capabilities and Methods of the X-12-ARIMA Seasonal Adjustment Program,”
Journal of Business and Economic Statistics, 16, 127-177.
[11] Fuller, W. (1987), Measurement Error Models, John-Wiley.
[12] Granger, C. W. J. and Hatanaka, M. (1964), Spectral Analysis of Economic
Time Series, Princeton U.P.
[13] Granger, C. W. J. and Newbold, P. (1977), Forecasting Economic Time Series,
Academic Press.
[14] Johansen, S. (1995), Likelihood Based Inference in Cointegrated Vector Autoregressive Models, Oxford UP.

28

[15] Harvey, A. and T. Trimbur (2008), ”Trend Estimation and The HodorickPrescott Filter,” J. Japan Statistical Society, 38-1, 41-49.
[16] Hayashi, F. (2000), Econometrics, Princeton University Press.
[17] Kitagawa, G. (2010), Introduction to Time Series Analysis, CRC Press.
[18] Kunitomo, N. , S. Sato and D. Kurisu (2018), Separating Information Maximum
Likelihood Estimation for High Frequency Financial Data, Springer.
[19] Kunitomo, N. and S. Sato (2017), ”Trend, Seasonality and Economic Time
Series : the Non-stationary Errors-in-variables Models,” SDS-4, MIMS, Meiji
University, http://www.mims.meiji.ac.jp/publications/2017-ds.
[20] Kunitomo, N., N. Awaya and D. Kurisu (2019), ”Some Properties of Estimation Methods for Structural Relationships in Non-stationary Errors-in-Variables
Models,” Japanese Journal of Statistics and Data Science, Springer (Open Access),[21] Morgenstern, Oskar (1950). On the accuracy of economic observations, (1966,
revised edition), Princeton University Press.
[22] Nerlove, M., D.M. Grether and J.L. Carvalho (1995), Analysis of Economic
Time Series : A Synthesis, Revised Edition, Academic Press.
[23] Müller, U. and M. Watson (2018), ”Long-run Covariability,” Econometrica,
86-3, 775-804.
[24] Nishimura, G. K., Seisho Sato, Akihiko Takahashi (2019), ”Term Structure
Models During the Global Financial Crisis: A Parsimonious Text Mining Approach,” Asic-Paciﬁc Financial Markets.
[25] Sato, S. and N. Kunitomo (2020a), ”On Backward Filtering for Noisy Nonstationary Time Series,” MIMS-RBP Statistics & Data Science Series (SDS-15),
MIMS, Meiji University.
[26] Sato, S. and N. Kunitomo (2020b), ”Frequency Regression and Smoothing for
Noisy Non-stationary Time Series,” In preparation.
[27] Yamada, H. (2019), ”A Smoothing Method that looks like the Hodrick-Prescott
Filter,” Unpublished Manuscript.

29

APPENDIX A : Mathematical Derivations
We now present some details of derivations that we have omitted in the previous
sections. We refer to Kunitomo, Sato and Kurisu (2018) as KSK(2018).
(i) On (4.24) and (4.27) : When we take m1 = 0 and m2 = m, then (4.27)
reduces (4.24). Hence, we show (4.27).
2π(j − 12 )(k − 12 ) (j, k = 1, · · · , n), we use the relation thatFor θjk = 2n+1
θjk + θj  ,k =

2π
1
2π
1


(j + j − 1)(k − ) , θjk − θj  ,k =(j − j )(k − ) .
2n + 12
2n + 1
2

We then have
(A.54)

4[cos θjk cos θj  ,k ]

(1)
k∈In



=

[e

i(θjk +θ  )
j ,k

+e

−i(θjk +θ 

j ,+k

)

]+

(1)
k∈In



[e

i(θjk −θ  )
j ,k

+e

−i(θjk −θ  )
j ,k

],

(1)
k∈In

(1)

(2)

where In = [1, · · · , m] (or In = [m1 + 1, · · · , m1 + m2 ]) is the index set for j andk.
(2)
For In = [m1 + 1, · · · , m1 + m2 ], by rewriting
1 2π
2π


θjk + θj  ,k = (m1 − )(j + j − 1) +(j + j − 1)(k − m1 ) ,2 2n + 1
2n + 1
and

1 2π
2π


(j − j ) +(j − j )(k − m1 ) ,θjk − θj  ,k = (m1 − )2 2n + 1
2n + 1the summation of the ﬁrst two terms in (A.54) becomese2πi(m1 + 12 ) 2n+1
(j+j −1)

2π

×1 − ei 2n+1 (j+j −1)m22π

1 − ei 2n+1 (j+j −1)

+e


2π
−i(m1 + 12 ) 2n+1
(j+j −1)

2π

×1 − e−i 2n+1 (j+j −1)m22π

1 − e−i 2n+1 (j+j −1)

For the last two terms in (A.54), we need to evaluate each term when (i) j = j and(2)(ii) j = j , separately. Using similar calculations in (A.54) with the index set In ,when j = j , the summation of last two terms becomese2πi(m1 + 12 ) 2n+1
(j−j )2π×

1 − ei 2n+1 (j−j )m22π



1 − ei 2n+1 (j−j )

+e



2π
−i(m1 + 12 ) 2n+1
(j−j )



2π

×



1 − e−i 2n+1 (j−j )m22π



1 − e−i 2n+1 (j−j ).

When j = j , θjk − θj  ,k = 0 and the summation of last two terms with the index set(2)
In becomes 2m2 . Hence, it is possible to evaluate each terms of (4.24) and (4.27).
30.

Using the relation
e

+e
2π

=

1



2π


2π
i(m1 + 12 ) 2n+1
(j+j −1)×

1 − ei 2n+1 (j+j −1)m22π

1 − ei 2n+1 (j+j −1)



2π



2π



2π
−i(m1 + 12 ) 2n+1
(j+j −1)×

1 − e−i 2n+1 (j+j −1)m2



1



2π

1 − e−i 2n+1 (j+j −1).

ei 2n+1 2 (j+j −1)(m1 ) − ei 2n+1 2 (j+j −1)(m1 +m2 )
2π



1

2π



1

ei 2n+1 (− 2 )(j+j −1) − ei 2n+1 ( 2 )(j+j −1)
2π

+



1



1

2π

e−i 2n+1 2 (j+j −1)(m1 ) − e−i 2n+1 2 (j+j −1)(m1 +m2 )
2π



1

2π

1



e−i 2n+1 (− 2 )(j+j −1) − e−i 2n+1 ( 2 )(j+j −1)

and the corresponding results for j − j (there are two cases when (a)j = j and(1)(b)j = j ), we obtain the results of (4.27), and then (4.24) by using In instead of(2)In .
(ii) The proof of Proposition 1 : Essentially, we apply the Central Limit
Theorem (Theorem 8.4.3 of Anderson (1971) or Theorem 7.6 of Durrett (1991)) to
the sequence of ergodic stationary (discrete) time series. We will give the basic stepsof our derivations and mention that the problem here is similar to those explainedin Chapters 7-9 of Anderson (1971) in details.
First, we need to show that the resulting variance-covariance terms correspondto those of the limiting Gaussian random variables.
For this purpose, we need to evaluate(A.55)
 
n


1

iθ
−iθ
(n)
(n) 
(eiθjk + e−iθjk )(e j k + e j k )E[rj rj  ] .
E zn (λk )zn (λk ) =2n + 1j,j =1When k = k , we ﬁnd that the right-hand side terms are bounded by using thestraightforward calculations. We notice that the right-hand side consists of sums offour terms associated withiθ(eiθjk + e−iθjk )(e j k + e(A.56)i(θjk +θ   )−iθ  j k−i(θjk +θ   ))i(θjk −θ   )

j k
j k
= e
+e
+e
= (A) + (B) + (C) + (D) (, say) .

j k

+e

−i(θjk −θ   )j kThen we ﬁnd that the sums of each terms associated with (A) and (B) in (A.55) arebounded, which becomes small when n is large. We write(A.57)nj,j =1ei(θjk +θ

jk

)



E[rj rj ] =

n−1




h=−(n−1) j  ∈S(h)

31

e

i(θ

 +θ   )h+j ,kj kΓ(h) ,

where S(h) = {1, 2, · · · , n − h} for h ≥ 0 and S(h) = {1 − h, 2 − h, · · · , n} for h < 0.
When h ≥ 0 given h, the second sum is given as
n−h




2π

1

1

j  =1

⎡
= ⎣



2π



1

1

ei 2n+1 (h+j − 2 )(k− 2 ) ei 2n+1 (j − 2 )(k − 2 ) × Γ(h)

n−h


⎤


2π



2π

1

1



1

ei 2n+1 (k+k −1)(j −1) ⎦ ei 2n+1 (h+1/2)(k− 2 )+ 2 (k − 2 ) × Γ(h) .



j =1

When h < 0 given h, the sum can be written
n+h−1




2π



1

1



1

1

ei 2n+1 [(h+j −1)(k+k −1)+ 2 (k− 2 )−(h+ 2 )(k − 2 )] × Γ(h) .

j  +h−1=0
Because these sums are ﬁnite and we have the condition +∞h=−∞ Γ(h) < +∞, thesum with (A) and (B) in (A.55) become arbitrarily small when n is large. (Theterms with (B) are the same as those with (A) except their signs in the exponential
parts.)
Second, we need to show that the sums of each terms with (C) and (D) in (A.55)
when k = k are dominant terms. We utilized the relation(A.58)nei(θjk −θ

jk

)

n−1




E[rj rj ] =





e

i(θ

 −θ   )h+j ,kj k

1× Γ(h) ,h=−(n−1) j  ∈S(h)j,j =1

When h ≥ 0 given h, the second sum is given as
n−h




2π



2π

1

1

1

ei 2n+1 (j −1)(k−k ) ei 2n+1 [(h+ 2 )(k− 2 )− 2 (k − 2 )] × Γ(h)j  =1and when h < 0 given h, the sum can be written
n+h−1


2π





1

1

1



1

ei 2n+1 [(h+j −1)(k−k )+ 2 (k− 2 )−(h− 2 )(k − 2 )] × Γ(h) .

j  +h−1=0The sums of the above terms are bounded when k = k by using the same argument as we have (A.57). (The terms with (D) are the same as those with (C)except their signs in the exponential parts.) On the other hand, when k = k ,n−h i2π k−k (j  −1)
e 2n+1= n − h. Then the dominant sums with (C) and (D) in (A.55)
j  =1
become n−1
k − 1/2nh [Γ(h) + Γ(−h)] .
cos 2π(A.59)
2n + 12n + 1h=−(n−1)32

We note that under the assumption of stationarity of rj , it has been known that
n
p

1


j =1 rh+j rj  −→ Γ(h) . (See Chapter 8 of Anderson (1971), and Brockwell andnDavis (1990), for instance.)(n)(n)
Third, we consider the situation that λk → s, λk → t as n → ∞ for 0 < s < t < 21 .
+∞Since h=−∞ Γ(h) < +∞ and Γ(h) is small as h is large, in the situation that(n)λk → λ for 0 < λ < 12 as n → ∞, we have (5.31). (We can take k such thatk/(2n) = λ + o(1/n ) ( > 0) such that |h|>n Γ(h) is arbitrary small.)
For the asymptotic normality, we set a sequence of random variables
 

(n)



Z n ek(n)(A.60)
Wk,k = (α1 , α2 ), k, k = 1, · · · , n,(n)
Z n ek and Zn = Pn Rn , Rn = (rj ) (n × p matrix), where αi (i = 1, 2) are p × 1 (non(n)zero) constant vectors and ek = (0, · · · , 0, 1, 0 · · · , 0) (k = 1, · · · , n) are n × 1 unit
vectors.
(n)
We can also construct a sequence of random variables Ws,t by using

(n)
(n)
ps,j = 2/(2n + 1) cos 2πs(j − 1/2) instead of pjk in Pn .
(n)

(n)

p

(n)

(n)

Then E[|Wk,k − Ws,t |2 ] → 0 as n → ∞ ( λk → s, λk → t (0 < s < t < 12 ).
Finally, we utilize the truncation argument used in the proof of Theorem 8.4.3
of Anderson (1971) and approximate (5.30) by the (ﬁnite order) moving-average
representation for a large h > 0
(n,h)
=μ+
ri

(A.61)hΓj wi−j

j=0

n
(n)
(n) (n,h)and zn,h (λk ) =. By applying the standard method of proof in
j=1 pjk rj
Theorem 8.4.3 of Anderson (1971), i.e., ﬁrst we ﬁx a h > 0 and we show the CLT(n)for zn,h (λk ) and then+∞by taking h → ∞ and we can show that the eﬀects of h are
negligible because s=−∞ Γ(s) < +∞. (Q.E.D.)
(x)

(x)

(s+v)

(iii) The proof of Proposition 2 : Let zk = (zkj ) and Zk
1, · · · , n) be the k-th row vector elements of n × p matrices(s+v)= (zkj) (k =

∗
(s+v)
Z(x)
= K∗n (Sn + Vn ) , K∗n = Pn C−1
n = Kn (Xn − Ȳ0 ) , Znn ,(A.62)

respectively, where we denote Xn = (xk ) = (xkg ), Sn = (sk ) = (skg ), Vn = (vk ) =

(x)
(s+v)
(vkg ), Zn = (zk ) (= (zkg )) are n × p matrices with zkg = zkg + zkg . We write
(x)

(s+v)

zkg , zkg , zkg

(x)

(s+v)

as the g−th component of zk , zk , zk

(k = 1, · · · , n; g = 1, · · · , p).

(f )
(x)
(x)
We use zkg (f = x, s+v) and decompose Σ̂x −Σx = (σ̂gh −σgh )gh ) for g, h = 1, · · · , p.

33

We re-writen1 (A.63)
zk zk − Σx
mn

 k=1
mn
mn

1 
1 
(x) (x)(s+v) (s+v)
=zk zk − Σx +E[zk zk]
mn k=1mn k=1
mn 
mn 


1 
1 
(s+v) (s+v)
(s+v) (s+v)
(x) (s+v)
(s+v) (x)zk z kzk zk.
− E[zk zk] ++ zk zk+
mn k=1mn k=1m

Then we will investigate the conditions that three terms except the ﬁrst one of(A.63) are op (1) andn1 p(x) (x)
zk zk − Σx −→ O
mn k=1m(A.64)

as mn → ∞ (n → ∞). Since Δxi is stationary, it is not diﬃcult to show the second
assertion 
via Chapter 5 of KSK(2018). By using (A.59), the expected (matrix) valuenof (1/m) m
k=1 zk zk is given by


n−1
n


k − 1/2n1
cos 2πh [Γ(h) + Γ(−h)] .
2n + 1m k=1
2n + 1h=−(n−1)

Then by Lemma 5.1 of KSK(2018), the middle sum is[1/2m] sin[2πmh/(2n + 1)]/ sin[πh/(2n + 1)],
which becomes 1 as m → ∞, m/n → 0+∞for a ﬁxed h. Then the limit is fSR (0) = −∞ Γ(h) because it is ﬁnite.
Here weuse the notation cij (i, j = 1, · · · , n) of Chapter 5 of KSK(2018) as cij =m(2/m) k=1 θik θjk . Then for (non-zero p × 1) constant vector α we can evaluate

E

mn
1 

mn k=1⎡2
(α zk )2 − α Σx α

= (2
)2 E ⎣
2n + 1≤ K1 [

2]2
2n + 1⎤2n

cjj  α (rj rj  − E(rj rj  ))α⎦j,j  =1nc2j,j  ,

j,j =1of fourth moments.
where K1 is a positive constant and we used the boundedness
nAgain by using Lemma 5.2 of KSK(2018), we ﬁnd that j,j  =1 c2j,j  = (n + 1/2)2 /m
and we have shown the ﬁrst assertion as m → ∞, m/n → 0.
Next we show that the second, third and fourth terms in the right-hand-side of(A.63) are asymptotically negligible as n → ∞. and we could estimate the variance34and covariance of the underlying processes consistently as if there were no seasonaland noise terms
(n)
= (0, · · · , 1, 0, · · · ) be an n × 1 vector.
Let bk = (bkj ) = ek Pn C−1
n = (bkj ) and ek
n(s+v)= j=1 bkj (sjg + vjg ) for the seasonal and noise part and use theWe write zkg
relation






−1
2
(Pn C−1n Cn Pn )k,k = δ(k, k )4 sin [(A.65)π1(k − )] .
2n + 12(s)(v)

Then under the conditions that Cj  = O(ρj ) (0 ≤ ρ < 1) and Cj  =O(ρj ) (0 ≤ ρ < 1), we can ﬁnd K2 (a positive constant) such that(s+v)
E[(zkg )]2 = E[

(A.66)n(s)bki (vig + vig )

i=1

n


(s)

bkj (vjg + vjg )] ≤ K2 × a∗kn .

j=1

It is because(s+v)
E[(zkg )]2 =

n


 (v)

(s)
bki bkj σgg(i − j) + σgg(i − j) ,

i,j=1
(v)

(s)

(s)

where σgg (i − j) + σgg (i − j) is the (i − j)−th autocovariance of vig + vig and(s)
vjg + vjg . (We denote bki = 0 for i < 0 and i > n.) Then(s+v)
E[(zkg )]2

=
≤

(s)

n−1


[

n


(v)
(s)
bkj bk,j+l (σgg(l) + σgg(l))]

l=−(n−1) j=1
n
∞


2
(v)
(s)
[
bkj ]][|σgg(l)| + |σgg(l)|] .
j=1
l=−∞
(v)

Because Cj  = O(ρj ) and Cj  = O(ρj ),Also it is straight-forward to ﬁnd that

∞

(v)
(s)
l=−∞ [|σgg (l)| + |σgg (l)|] is bounded.

mnmn
2k − 11 
1 m2∗
1 − cos(πakn =2) = O( 2n ) ,mn k=1mn k=12n + 1n
by using the relation
m


m
2π

sin( 2n+1m)2π2π2k − 1i 2n+1(k− 12 )−i 2n+1(k− 12 )2 cos(π
[e
+e
]=
)=
.
π
2n + 1sin( 2n+1 )k=1
k=1

Then the second term of (A.63) becomes(A.67)

mn
mn
m2
1 
1 
(s+v) 2E[zkg ] ≤ K3
a∗kn = O( 2n ) ,mn k=1mn k=1n

35

which is o(1) if we set α such that 0 < α < 1 and K3 is a positive constant.
(The arguments here are similar to the derivations in Chapter 5 of KSK(2018), butthere is a major diﬀerence on the conditions because their asymptotic arguments
are based on the situation when
√ the length of observed intervals decrease. In thepresent situation, there is no n factor in Kn -transformation and we use K∗n anda∗kn while they have used akn = na∗kn .)For the fourth term,
2

mn
mn


1 
mn
1 
(x) (v)(x) (x) (v) (v)zkg zkg= 2
E zkg zk ,g zkg zk ,g = O( 2 ) .
E
mn j=1mn nk,k =1

2π
In the above evaluation we have used the evaluation that if we set sjk = cos[ 2n+1
(j −1
1)(k − 2 )] (j, k = 1, 2, · · · , n), then we have the relation
2

|

n


sjk sj,k | ≤ [

j=1

n


s2jk ] =j=1

n 1+ for any k ≥ 1 .
2 4(See Chapter 5 of Kunitomo, Sato and Kurisu (2018).)
Finally, for the third term, we need to consider the variance of(s+v) 2

=

(zkg
n


j,j  =1

(s+v)

) − E[(zkg )2 ]


(s)
(s)
(s)
(s)
bkj bk,j  (vjg + vjg )(vj  ,g + vj  ,g ) − E[(vjg + vjg )(vj  ,g + vj  ,g )] .

Then by using the assumption we have after some evaluations, we ﬁnd a positive
constant K4 such that
2

mn
1 
(s+v) 2(s+v)
(A.68) E((z) − E[(zkg )2 ])mn k=1 kg

mn
n


1 
(s)
(s)
=
E
b
b
k1 ,j1 k1 ,j2 (vj1 ,g + vj1 ,g )(vj2 ,g + vj2 ,g )m2n k ,k =1
j1 ,j2 ,j3 ,j4 =1
1 2

(s)
(s)
−E(vj1 ,g + vj1 ,g )(vj2 ,g + vj2 ,g )


(s)
(s)
(s)
(s)
×bk2 ,j3 bk2 ,j4 (vj3 ,g + vj3 ,g )(vj4 ,g + vj4 ,g ) − E(vj3 ,g + vj3 ,g )(vj4 ,g + vj4 ,g )
mn
1
m3
1 
a∗kn ]2 = O( 2 × ( 2n )2 ) ,≤ K4 2 [mn k=1mnnwhich is O(m4n /n4 ) by straight-forward calculations.
Here we just give an illustration of our derivations when p = 1 and si = 0 (i =36

1, · · · , n). We re-write (3.12) as vi =((v)
j=0 cj ei−j and we evaluate

m

1 2 
bk ,j bk ,j bk ,j bk ,j)m k ,k =1 j ,j ,j ,j 1 1 1 2 2 3 2 4
1

= (

∞

2

1

2

3

4

×E{[vj1 vj2 − E(vj1 vj2 )][vj3 vj4 − E(vj3 vj4 )]}∞
bk1 ,j1 bk1 ,j2 bk2 ,j3 bk2 ,j4
c l1 c l 2 c l 3 c l 4

m


1 2
)
m k ,k =1 j ,j ,j ,j
1

2

1

2

3

l1 ,l2 ,l3 ,l4 =04×E{[ej1 −l1 ej2 −l2 − E(ej1 −l1 ej2 −l2 )][ej3 −l3 ej4 −l4 − E(ej3 −l3 ej4 −l4 )]} .
We need to evaluate the corresponding terms for four cases when (i) j1 −l1 = j2 −l2 =j3 −l3 = j4 −l4 , (ii) j1 −l1 = j2 −l2 = j3 −l3 = j4 −l4 , (iii) j1 −l1 = j3 −l3 = j2 −l2 =j4 − l4 , (iv) j1 − l1 = j4 − l4 = j2 − l2 = j4 − l4 . By using the condition in the general

(s)
(v)
case that Cj  = O(ρj ), Cj  = O(ρj ) (j ≥ 0, 0 ≤ ρ < 1), we have ∞j=0 |cj | < ∞n
in this special case. We also utilize the relation such as j=1 bkj bk j = δ(k, k )a∗kn
in the general case and we have used the notation that bk,j = 0 for k = 1, · · · , m,
j < 0, j > n and cj = 0 (j < 0).
Then in each (i)-(iv) case, we can take a positive constant K5 such that (A.68) is
less than
n
n
n
n
m



1 2   2
1/2
2
1/2
2
1/2
K5 ( )[
bk1 ,j1 ] [bk1 ,j2 ] [bk2 ,j3 ] [b2k2 ,j4 ]1/2 .
m k ,k =1 j =1
j =1
j =1
j =1
1

2

1

2

3

4

The above evaluation methd works in the general case with a complication of notations. Therefore, by using (68), the third term of (A.63) is negligible if we set α
such that 0 < α < 1. (Q.E.D)

APPENDIX B : Some Figures
We now present ﬁgures used in Section 7. As we have explained in Section 2, alldata are oﬃcial data published by ESRI (Economic and Social Research Institute),
Cabinet Oﬃce of Japan, METI, and Statistics Bureau, Ministry of Internal Aﬀairs
and Communications. They are available from the government oﬃcial website :e.stat.

37

6e+05
4e+05
2e+05

RCON

0e+00
−2e+05
0.0

0.1

0.2

0.3

0.4

0.5

freq

Figure 5: Orthogonal Process zk for Original Series. (Without diﬀerencing C−1n

RCON−10000−5000

0

we used the transformation Pn of data and then calculated zk (k = 1, · · · , n). Data are thequarterly real consumption between 1994Q1-2018Q2, published in 2018 by the Economic Social
Research Institute (ESRI), Cabinet Oﬃce, Japan.)

0.0

0.1

0.2

0.3

0.4

0.5

freq

Figure 6: Orthogonal Process zk for Quarterly Consumption.
(We used
the transformation Kn of data and then calculated zk (k = 1, · · · , n). Data are the quarterly
real consumption between 1994Q1-2018Q2, published in 2018 by the Economic Social ResearchInstitute (ESRI), Cabinet Oﬃce, Japan.)

38

1.0
0.8
0.6

RCON

0.4
0.2
0.0
0.0

0.1

0.2

0.3

0.4

0.5

freq

Figure 7: Squared Orthogonal Process z2k for Quarterly Consumption. (We

2000
−2000

0

RCON

4000

6000

calculated the cumulative sum of squares of the transformed orthogonal z2n from the data.)

0.0

0.1

0.2

0.3

0.4

0.5

freq

Figure 8: Orthogonal Process zk for Seasonally Diﬀerenced Consumption.

(After seasonal diﬀerencing instead of Cn−1 , we used the transformation Pn on the diﬀerenced dataand calculated the orthogonal elements.)

39

Figure 9: Seasonally Ajusted Series. (The oﬃcial seasonally adjusted (S.A.) series areconstructed using X-12-ARIMA and we compared the SIML ﬁltered value. Data are the quarterly
real GDP between 1994Q1-2018Q2, published in 2018 by the Economic Social Research Institute(ESRI), Cabinet Oﬃce, Japan.)Figure 10: Monthly Consumption Series. (We compared the original series and the
SIML ﬁltered value. Data are the monthly real consumption between 2002M1-2016M12, publishedin 2017 by the Statistics Bureau, Japan.)40

Figure 11: Monthly Consumption Series. (We compared the original series and the
SIML ﬁltered value. Data are the monthly consumptions between 2002M1-2016M12, published in2017 by Ministry of Economy, Trade and Industry (METI), Japan.)Figure 12: Consumption Series. (Data are the quarterly real consumption between
1994Q1-2018Q2, published in 2018 by the Economic Social Research Institute (ESRI), Cabinet
Oﬃce, Japan.)

41

Figure 13: Predicted Series. (The predicted and monthly consumption state variable areshown between 1994Q1-2018Q2.)42