科研費シンポジウム「統計科学と関連分野における諸問題に関する理論と方法論の革新的展開」（Innovative developments in theory and methodology regarding various issues in statisticalscience and related ﬁelds）
日時：2023 年 9 月 29 日（金）〜10 月 1 日（日）(Date: Friday, 29, September – Sunday, 1, October, 2023)
場所：新潟大学駅南キャンパスときめいと 講義室 A,B(Place: Niigata University, Station South Campus Tokimate, Lecture Room A, B)（TEL: 025-248-8141）
科学研究費・基盤研究（A) （課題番号：20H00576）「大規模複雑データの理論と方法論の革新的展開」
研究代表者：青嶋 誠（筑波大学）

Makoto Aoshima (University of Tsukuba)
開催責任者：蛭川 潤一（新潟大学）

Junichi Hirukawa (Niigata Universiry)Program
Friday, 29, September
Reception 13:00-13:15
Opening 13:15-13:20 Junichi HIRUKAWA (Niigata University)
Afternoon Session I (in Japanese) 13:20-16:00
Chair Shimatani, K. Ichiro (The Institute of Statistical Mathematics)
1. 13:20-14:00 江頭 健斗, 矢田 和善, 青嶋 誠
東京理科大学創域理工学部, 筑波大学数理物質系, 筑波大学数理物質系

Hierarchical clustering and its asymptotic behaviors in multiclass HDLSS settings
2. 14:00-14:40 齊藤 実祥, 寒河江 雅彦
金沢大学 人間社会研究域 経済学経営学系
多次元 Tensor Product 型密度関数の推定とその漸近的性質

3. 14:40-15:20 種市 信裕, 関谷 祐里
元北海道教育大札幌校, 北海道教育大学釧路校

2nd order correction term を最小化する分割表のΦ‐ダイバージェンス独立性検定統計量
4. 15:20-16:00 柿沢 佳秀
北海道大学大学院経済学研究院
境界バイアスを回避する非対称カーネル法の様々な適用について

Coﬀee Break 16:00-16:30
Afternoon Session II (in Japanese & English) 16:30-17:50
Chair Yuichi Goto (Kyushu University)
5. 16:30-17:10 曾 小強 (Xiaoqiang Zeng), 柿沢 佳秀 (Yoshihide Kakizawa)
Hokkaido University
Some estimators in the ADCINAR(1) process
6. 17:10-17:50 竹内 努 (Tsutomu T. TAKEUCHI)
名古屋大学素粒子宇宙物理学専攻 (Division of Particle and Astrophysical Science, Nagoya Uni-

versity, Japan)
多様体学習による銀河進化の定量化と定式化 (Quantiﬁcation and Formulation of Galaxy Evolution

by Manifold Learning)Saturday, 30, September
Morning Session (in Japanese) 9:30-10:50
Chair Yan Liu (Waseda University)
7. 9:30-10:10 永井 勇
中京大学 教養教育研究院

GMANOVA モデルにおけるフルランク仮定が不要な罰則無推定法
8. 10:10-10:50 藤森 洸, 佃 康司
信州大学, 九州大学

Two step estimations via the Dantzig selector for models of stochastic processes with highdimensional parameters
Lunch 10:50-13:20
Afternoon Session I (in English) 13:20-15:20
Chair Kou Fujimori (Shinshu University)
9. 13:20-14:00 島谷健一郎 (Shimatani, K. Ichiro)
統計数理研究所 (The Institute of Statistical Mathematics)
空間種分布モデルに対する多様性指数の統計数理 (Statistical mathematics of biodiversity indicesfor spatial species distribution model)
10. 14:00-14:40 張 元宗 (Chang Yuan-Tsung), 篠崎 信雄 (Nobuo Shinozaki)
目白大学 (Mejiro university), 慶應大学 (Keio university)

Simultaneous estimation of Poisson means: recent developments and its applications
11. 14:40-15:20 田中 勝人 (Katsuto Tanaka)
Hitotsubashi University (Emeritus Professor)
Extensions of Darling’s Formula for the Fredholm DeterminantCoﬀee Break 15:20-15:50

Afternoon Session II (Guest speakers session) 15:50-17:50
Chair Junichi Hirukawa (Niigata University)
12. 15:50-16:30 Liang-Ching Lin
National Cheng Kung University
Authors: Liang-Ching Lin (National Cheng Kung Univ.), Hsiang-Lin Chien (National Cheng
Kung Univ.), Hao Sung (National Cheng Kung Univ.) and Sangyeol Lee (Seoul National University)
Comprehensive interval-valued time series model with application to the S&P 500 index andPM2.5 level data analysis
13. 16:30-17:10 Nan-Jung Hsu
National Tsing Hua University
Authors: Nan-Jung Hsu (National Tsing-Hua Univ), Hsin-Cheng Huang (Institute of Statistical
Science, Academia Sinica), Ruey S. Tsay ( University of Chicago) and Tzu-Chieh Kao (National
Tsing-Hua Univ.)
Matrix Autoregressive Spatio-Temporal Models
14. 17:10-17:50 Hsin-Cheng Huang
Institute of Statistical Science, Academia Sinica
Authors: Chien-Chung Wang (Colorado State Univ.), Chih-Hao Chang (National Chengchi
Univ.) and Hsin-Cheng Huang (Institute of Statistical Science, Academia Sinica)
Stable p-value Assignment in High-Dimensional Regression via Data Splitting
Conference Dinner 18:30-21:30
いかの墨 新潟駅前店(Ika-no-sumi)
6.000 yen
(TEL: 025-242-0510)

Sunday, 1, October
Morning Session (in English) 9:30-12:10
Chair Katsuto Tanaka (Hitotsubashi University)
15. 9:30-10:10 塩濱 敬之 (Takayuki Shiohama)
南山大学理工学部 (Department of Data Science, Nanzan University)

Modeling Joint Cylindrical Distributions and Related Markov Processes
16. 10:10-10:50 後藤 佑一 (Yuichi Goto)
九州大学 (Kyushu University)

Authors: Yuichi Goto (Kyushu Univ.), Xuze Zhang (Maryland Univ.), Benjamin Kedem (Maryland Univ.) and Shuo Chen (Maryland Univ.)
Test for the existence of the residual spectrum with application to brain functional connectivity
detection
17. 10:50-11:30 Yan Liu
Faculty of Science and Engineering, Waseda University
Prediction-based statistical inference for multiple time series
18. 11:30-12:10 Junichi Hirukawa, Kou Fujimori
Niigata University, Shinshu University
Weak convergence of the partial sum of I(d) process to a fractional Brownian motion in ﬁnite
interval representation
Closing 12:10-12:15 Junichi HIRUKAWA (Niigata University)

Hierarchical clustering and its asymptotic behaviors
in multiclass HDLSS settings
Kento Egashiraa , Kazuyoshi Yatab , Makoto AoshimabaFaculity of Science and Technology, Tokyo University of SciencebInstitute of Mathematics, University of Tsukuba

Hierarchical clustering is a methodology to group a set of datas by building dendrogram based on a similarity or a dissimilarity between clusters so that datas in acluster are similar in the sense of pre-determined linkage function. In hierarchicalclustering, one can observe a process how a cluster is combined or devided throughdendrogram on graphic. Hierarchical clustering has been approved as useful tool
for analysis of gene expression microarray data. In fact, applications of hierarchicalclustering on gene expression microarray data are given by Eisen et al. [4], Perouet al. [8], Bhattacharjee et al. [2], among others. A characteristic of datas used in
Eisen et al. [4], Perou et al. [8] and Bhattacharjee et al. [2] is that the number ofvariables is much larger than sample size. This type of data represented by gene expression microarray data is called high-dimension, low-sample-size (HDLSS) data.
Substantial work about clustering has been performed on HDLSS asymptotics inrecent years. For example, Liu et al. [6] proposed a two-way split clustering called“statistical signiﬁcance of clustering(SigClust)” especially for HDLSS data. Ahn
et al. [1] proposed a hierarchical divisive clustering and considered its high dimensional asymptotics. Huang et al. [5] developed the SigClust by Liu et al. [6] with soft
thresholding approach. Yata and Aoshima [9] gave consistency properties of sample principal component scores and applied it to clustering under high dimensionalsettings. Nakayama et al. [7] investigated clustering by kernel principal component
analysis for HDLSS data. Borysov et al. [3] studied behaviors of hierarchical clustering under several asymptotic settings from moderate dimension through HDLSS,
nevertheless it is considered that theoretical assumptions are strict for HDLSS data
due to having discussions on several asymptotic settings at once.
Given this background, we focused on HDLSS settings and considered asymptotic properties of hierarchical clustering with several linkage functions. This study
explores practical assumptions to indicate the behavior of hierarchical clustering.
In this talk, we theoretically investigated hierarchical clustering when both thedimension and sample size approach inﬁnity at ﬁrst. Then, we gave an asymptotic
behavior in boundary cases of thresholds which divide asymptotic behaviors. Finally,
we examined the hierarchical clustering theoretically in multiclass HDLSS context
and reported numerical simulations.

References[1] Ahn, J., Lee, M.H., Yoon, Y.J. (2012). Clustering high dimension, low sample
size data using the maximal data piling distance. Statistica Sinica, 22, 443–464.
[2] Bhattacharjee, A., Richards, W. G., Staunton, J., Li, C., Monti, S., Vasa, P.,
Ladd, C., Beheshti, J., Bueno, R., Gillette, M., Loda, M., Weber, G., Mark,E. J., Lander, E. S., Wong, W., Johnson, B. E., Golub, T. R., Sugarbaker, D.
J., Meyerson, M. (2001). Classiﬁcation of human lung carcinomas by mRNA
expression proﬁling reveals distinct adenocarcinoma subclasses. Proceedings of
the National Academy of Sciences of the United States of America, 98, 13790–
13795.
[3] Borysov, P., Hannig, J., Marron, J.S. (2014). Asymptotics of hierarchical clustering for growing dimension. Journal of Multivariate Analysis, 124, 465–479.
[4] Eisen, M. B., Spellman, P. T., Brown, P. O., Botstein, D. (1998). Cluster
analysis and display of genome-wide expression patterns. Proceedings of the
National Academy of Sciences of the United States of America, 95, 14863–
14868.
[5] Huang, H., Liu, Y., Yuan, M., Marron, J.S. (2015). Statistical Signiﬁcanceof Clustering using Soft Thresholding. Journal of computational and graphicalstatistics, 24, 975–993.
[6] Liu, Y., Hayes, D.N., Nobel, A., Marron, J.S. (2008). Statistical signiﬁcance ofclustering for high-dimension, low-sample size data. Journal of the American
Statistical Association, 103, 1281–1293.
[7] Nakayama, Y., Yata, K., Aoshima, M. (2021). Clustering by principal component analysis with Gaussian kernel in high-dimension, low-sample-size settings.
Journal of Multivariate Analysis, 185, 104779.
[8] Perou, C. M., Sørlie, T., Eisen, M. B., van de Rijn, M., Jeﬀrey, S. S., Rees,
C. A., Pollack, J. R., Ross, D. T., Johnsen, H., Akslen, L. A., Fluge, O.,
Pergamenschikov, A., Williams, C., Zhu, S. X., Lønning, P. E., Børresen-Dale,
A. L., Brown, P. O., Botstein, D. (2000). Molecular portraits of human breast
tumours. Nature, 406, 747–752.
[9] Yata, K., Aoshima, M. (2020). Geometric consistency of principal componentscores for high-dimensional mixture models and its application. Scandinavian
Journal of Statistics, 47, 899–921.

多次元 Tensor Product 型密度関数の推定とその漸近的性質
齊藤実祥 (金沢大学), 寒河江雅彦 (金沢大学)

1 研究目的多次元ヒストグラムにおいて、各ビン間の不連続性によって推定精度が低下する問題がある。そこで、平滑化によるモデル拡張でヒストグラムの推定精度 O(n−2q/2q+d ), (q = 1), (q は推定効率、d は次元数) の改良を考える。先行研究では、Scott(1985) と Hjort(1986) がヒストグラムの各ビンの中点を線形で接続する

Frequency Polygon(FP) を d 次元へ拡張し、その推定精度が O(n−4/d+4 ), (q = 2) で、ヒストグラムを改良することを示した。本稿では、曲面の近似で用いられる Tensor Product Spline を密度関数の推定に適用し、その理論面を整備する。2 次元の場合での Tensor Product Spline 密度推定量 (以降、TPS 推定量) について、
平均積分二乗誤差 (MISE) に基づく推定精度と漸近正規性を示す。その中で、TPS 推定量で推定精度の維持とパラメータの削減が両立できる方法について提案する。また、TPS 推定量の改良モデルを提案し、その漸近的性質から推定精度が改良できることを示す。

2 Tensor Product Spline 密度推定量
2.1 TPS 推定量の構築と漸近的性質
R2 上の確率密度関数 f (x, y) から n 個のデータを与え、ヒストグラムのビンの中点を x 軸側 xi , (i =1, . . . , N )、y 軸側 yj , (j = 1, . . . , M )、ビン幅を x 軸側 h1 、y 軸側 h2 、i, j 番目のビン Bi,j の定義域を[] []Bi,j = xi − h21 , xi + h21 × yj − h22 , yj + h22 とする。
ビン Bi,j において、提案モデル Tensor Product Spline 密度推定量 (以降、TPS 推定量)fˆi,j (x, y) は次の通
りである;

(

fˆi,j (x, y) = 1 (x − xi )
[

(i,j)

ただし、A = akl

]


a11
)
(x − xi )2 a21
a31

a12
a22
a32



a13
1
a23   (y − yj )  ,a33(y − yj )2(x, y) ∈ Bi,j ,(1), (k, l = 1, 2, 3) は未知の定数で、簡単のため各要素を akl と略す。

未知の定数 akl を求める方程式について次の行列表現が可能で、一意に導出できる；

0

A =  12− 2h1 1

1
0
0

0



γi−1,j−1
1 βi,j−12
1γi,j−12h1

αi−1,j
ui,j
αi,j


0
γi−1,j


βi,j
1
γi,j
0

1
2

0
1
2


− 2h1 2
0 ,

(2)

1
2h2

ただし、既知の各係数について、Bi,j の中点での高さは度数 νi,j を用いて ui,j = νi,j /nh1 h2 、x での偏微分˙˙¨α·,· = fˆx (·, ·)、y での偏微分 β·,· = fˆy (·, ·)、x, y 両方での偏微分 γ·,· = fˆxy (·, ·) である。
TPS 推定量 fˆ(x, y) の漸近的な MISE(AMISE) は、{} 1.1125h4 ( )h4 ( ) h2 h2 ( )AMISE fˆ(x, y) =+ 1 R f¨xx + 2 R f¨yy + 1 2 R f¨xy ,(3)nh1 h2
576
576
288
∫
ただし、R(ϕ) = ϕ(x, y)2 dxdy である。
()()
このとき、最適ビン幅は h∗i = O n−1/6 、最小 AMISE は O n−2/3 , (q = 2) で 2 次元 FP と漸近同等である。また、fˆ(x, y) は漸近正規性が成り立つ。

2.2 TPS 推定量のパラメータ削減法
ここで、A でパラメータを全 9 個から 6 個に減らす改良法を提案する。この理由は、多項式における 3 次
以上の項 (x − xi )k−1 (y − yj )l−1 , k + l ≥ 5 に対応する 3 項を消去することで、MISE の収束レートは変えずに定数項を改良できるからである。

パラメータを削減した TPS 推定量は次の通り表現される；a11)(x − xi )2 a21
a31(

fˆi,j (x, y) = 1 (x − xi )

a12
a22
0



a13
1
0   (y − yj )  ,0(y − yj )2(x, y) ∈ Bi,j .

(4)

この時の AMISE は次の通りである；{} 1.0696
h4 ( )h4 ( ) h2 h2 ( )AMISE fˆ(x, y) =+ 1 R f¨xx + 2 R f¨yy + 1 2 R f¨xy .
nh1 h2
576
576
288

(5)

また、パラメータを削減した fˆ(x, y) についても漸近正規性が成り立つ。

3 TPS 推定量の改良
各ビンで面積相等性を満たすように制約条件を追加することで TPS 推定量を改良することを考える。
ビン Bi,j における改良版 TPS 推定量 ĝi,j (x, y) は次の多項式で表現され、制約条件に基づく方程式を解く
ことで各係数が導出される；

ĝi,j (x, y) = a31 (x − xi )2 + a13 (y − yj )2 + a22 (x − xi )(y − yj )+ a21 (x − xi ) + a12 (y − yj ) + a11 .

(6)

改良版 TPS 推定量 ĝ(x, y) の漸近的な MISE(AMISE) は、(... ))
4.02101 ( 6 (... )h R f xxx + h62 R f yyy+
nh1 h230240 1h4 h2 (... ) h2 h4 (... )(7)+ 1 2 R f xxy + 1 2 R f xyy .
6912
6912
(
)
(
)
このとき、最適ビン幅は h∗i = O n−1/8 、最小 AMISE は O n−3/4 , (q = 3) であり、TPS 推定量の推AMISE(ĝ(x, y)) =

定精度を改良することが示された。
また、改良版 TPS 推定量 ĝ(x, y) は漸近正規性が成り立ち、h ∝ O(n−α ), (x, y) ∈ Bi,j に対して、√
1dのとき： nh1 h2 {ĝi,j (x, y) − g(x, y)} −→ N (Bias[ĝi,j (x, y)], 4.0210g(ζi , ηj )) ,8√1dα > のとき： nh1 h2 {ĝi,j (x, y) − g(x, y)} −→ N (o(1), 4.0210g(ζi , ηj )) ,8∫
ただし、g(ζi , ηj ) は pi,j = Bi,j g(u, v)dudv = h1 h2 g(ζi , ηj ), (ζi , ηj ) ∈ Bi,j を満たす点である。
α=

(8)
(9)

4 結論
多次元ヒストグラムを Tensor Product 型で平滑化する 2 次元 TPS 推定量の推定精度は AMISE(fˆ(x, y)) =()()
O n−2/3 , (q = 2) で、2 次元ヒストグラムの O n−1/2 , (q = 1) を改良できるものの、2 次元 FP の()
O n−2/3 , (q = 2) と同等の収束レートである。2 次元 TPS 推定量でパラメータを全 9 個から 6 個に削減す()
る方法を提案し、収束レートが O n−2/3 , (q = 2) で変化せず、分散の定数項が低下することで AMISE が改良できることを示した。また、面積相等性の条件を追加して TPS 推定量を改良することを提案した。その()

推定精度が AMISE(ĝ(x, y)) = O n−3/4 , (q = 3) で、TPS 推定量より優れた推定精度であることを示した。

参考文献[1] Hjort, N. L. (1986), “On Frequency Polygons and Average Shifted Histograms in Higher Dimensions”, Tech Report 22, Stanford University.
[2] Scott, D. W. (1985), “Frequency Polygons: Theory and Application”, Journal of the American
Statistical Association, Vol. 80, No. 390, pp.348-354.

2nd order correction term を最小化する分割表の ϕ-ダイバー
ジェンス独立性検定統計量
元北海道教育大札幌校 種市信裕
北海道教育大釧路校 関谷祐里

1. ϕ-ダイバージェンスに基づく分割表独立性検定統計量．
多次元の J1 × J2 × · · · × JM 分割表において，サンプルサイズ n が固定されている
多項分布モデルの場合を考える。すなわち，(j1 , . . . , jM ) セル (jm = 1, . . . , Jm ; m =1, . . . , M ) の観測度数を Xj1 ,...,jM とし，X ∗ = (X11···1 , X11···12 , . . . , X11···1JM , . . . ,∑∑ M∗
XJ1 J2 ···JM )′ とおく。ただし， Jj11=1 · · · JjM=1 Xj1 ···jM = n である。すると，X は多項分布 MultK (n, p∗ ) に従って分布する。ただし，K は分割表のセル数，すなわ∏∗′
ちK = M
m=1 Jm , であり，p = (p11···1 , p11···12 , . . . , p11···1JM , . . . , pJ1 J2 ···JM ) , 0 <∑J1∑JM
pj1 ···jM < 1, (jm = 1, . . . , Jm ; m = 1, . . . , M ), j1 =1 · · · jM =1 pj1 ···jM = 1 で(I)ある。M 次元分割表の完全独立性の仮説は，H0 : pj1 ···jM = P (j1 , . . . , jM ), (jm =∏1, . . . , Jm ; m = 1, . . . , M ) である。ただし，P (j1 , . . . , jM ) = M
m=1 p·(m,jm ) , (jm =1, . . . , Jm ; m = 1, . . . , M ), であり，記号 a·(m,jm ) は aj1 ···jm ···jM , (jm = 1, . . . ,∑Jm ; m = 1, . . . , M ) を，M 個の添え字をもつ列とするとき，a·(m,jm ) = Jj11=1∑JM
jm+1 =1 · · ·jM =1 aj1 ···jm ···jM , (jm = 1, . . . , Jm ; m = 1, . . . , M )(I)で与えられる。H0 を検定するための ϕ-divergence に基づく統計量 Cϕ は，Cϕ =∑∑ M
2n Jj11=1 · · · JjM=1 P̂ (j1 , . . . , jM )ϕ(p̂j1 ···jM /P̂ (j1 , . . . , jM )) で定義される。ここで，ϕ(t)は t > 0 で定義され ϕ(1) = ϕ′ (1) = 0, ϕ′′ (1) = 1 (Pardo et al. [2]) を満たす実凸関

···

∑Jm−1

∑Jm+1

jm−1 =1∏数であり，P̂ (j1 , . . . , jM ) = M
m=1 p̂·(m,jm ) , p̂j1 ···jM = Xj1 ···jM /n, (jm = 1, . . . , Jm ;m = 1, . . . , M ), p̂·(m,jm ) = X·(m,jm ) /n, (jm = 1, . . . , Jm ; m = 1, . . . , M ) である。
(I)

Zografos [4] は M = 2 のときにこの統計量を示した。完全独立性の仮説 H0 のもと∑で，統計量 Cϕ はすべて，自由度 η = K − M
m=1 Jm + M − 1 の極限カイ二乗分布をもつ。ϕ(t) として，関数 ϕa (t) = {a(a + 1)}−1 {ta+1 − t + a(1 − t)}, (a ̸= 0, −1); =t ln t + 1 − t, (a = 0); = − ln t − 1 + t, (a = −1) を選ぶと，Cϕa は，Cressie & Reada になる。R0 は対[1] により提案されたパワーダイバージェンスに基づく統計量 R(I)(I)
2/3

1 は Pearson X 2 統計量，R数尤度比統計量，R(I)(I) は Cressie & Read [1] によってQ推奨された統計量である。ϕ として ϕa (t) = (t − 1)2 /[2 {a + t(1 − a)}], a ∈ [0, 1] を選ぶと Rukhin [3] により提案された統計量 Qa(I) となる。Q1(I) は Pearson X 2 統計量である。

2. 2nd-order correction term．
(I)

統計量 Cϕ の帰無仮説 H0 のもとでの原点のまわりの s 次モーメントの評価のため(I)−1の展開式 E((Cϕ )s |H0 ) = E((Fη )s ) + mC
(s = 1, 2, . . .) を考える。
ϕ (s)/n + o(n ),ここで，Fη は，自由度 η のカイ二乗分布に従う確率変数とする。この時，mC
ϕ (s) は，C
2nd-order correction term と呼ばれる。もし，mϕ (s) が０に近ければ，その統計量の分布が自由度 η のカイ二乗分布に近いと考えることができる。故に，mC
ϕ (s) の値を計算することによって，どの統計量がカイ二乗分布に近いかを調べることができる。M 次元分割表における完全独立性検定統計量 Cϕ はモーメントが存在しない場合がある。しかしながら，(Cϕ )s = ∞ になる確率は n の増加に対して非常に速く０に近づく。つまり，大きなサンプルサイズ n に対して，Cϕ の n−1 のオーダーまでの展開
式を Cϕ∗ とおくとき，Cϕ の分布は，Cϕ∗ の分布に非常によく似ていると考えられる。
よってここでは，s 次モーメントの評価のための展開式において Cϕ のかわりに Cϕ∗ とおいた場合の 2nd-order correction term を mC
ϕ (s) とする。

3. 多次元分割表における完全独立性検定の場合の 2nd-order correction
term に関する定理および導出される検定統計量.
mC
ϕ に関して次の定理を得る。
定理 1. Cϕ∗ の分布は連続分布で，ϕ(t) を４回微分可能，ϕ(4) (t) が t = 1 で連続で，
pj1 ···jM = O(K −1 ), (jm = 1, . . . , Jm ; m = 1, . . . , M ) とする。このとき，ϕ-ダイバージェンス統計量に対する展開式 Cϕ∗ に対して，mC
ϕ (s) = 0, (s = 1, 2, . . .) を満たす ϕ(I)について，H0 のもとで，J1 , . . . , JM が同程度のオーダーでかつ K =無限大に発散するとき，4ϕ′′′ (1) + 3ϕ(4) (1) は０に収束する。

∏M

m=1 Jm がmC
ϕ (s) の評価に関して以下の定理を得る。
定理 2. Cϕ∗ の分布は連続分布で，ϕ(t) を４回微分可能，ϕ(4) (t) が t = 1 で連続で，pj1 ···jM = O(K −1 ), (jm = 1, . . . , Jm ; m = 1, . . . , M ) とする。このとき，4ϕ′′′ (1) +
C
′′′
s ∑M
3ϕ(4) (1) = 0 を満たす ϕ に対して，mC
m=1 Sm /Jm
ϕ (s) は，mϕ (s) = −(s/3)ϕ (1)K
∑M
2
s
s−1
+O(K ) +O(K ( m=1 Jm ) ), (s = 1, 2, . . .) と評価される。
定理 1, 定理 2 より，より速くカイ二乗分布に収束する検定統計量を得るために，
4ϕ′′′ (1) + 3ϕ(4) (1) = 0 かつ ϕ′′′ (1) = 0，つまり，ϕ′′′ (1) = ϕ(4) (1) = 0 を満たす関数
ϕ がもし存在するならば，その ϕ により定義される検定統計量 Cϕ を推奨する。このような性質は，パワーダイバージェンス統計量および Rukhin による統計量の場合に1 , Q1 (共に Pearson X 2 統計量）のみが満たす。一方，ψ (t) = (3/2)ϕ (t) −は，R(I)aa(I)Q(1/2)ϕa (t) とおき，a = 2/3 の場合を ϕN T (t) とする，つまり，ϕN T (t) = ψ2/3 (t) =(27/20)t5/3 − 3t − (27/4)(t + 2)−1 + 39/10 とおき，これを ϕ とすると上述の性質，つまり「ϕ(1) = ϕ′ (1) = ϕ′′′ (1) = ϕ(4) (1) = 0, ϕ′′ (1) = 1, かつ t > 0 において凸」を満たす。よってこの ϕN T を用いることによって，推奨される検定統計量を構築することができる。本報告では，構築された統計量と Pearson X 2 統計量を含めた他の検定統計量の性能を検討する。

参考文献[1] Cressie, N., and Read, T. R. C. JRSS, B, 46, (1984), 440–464.
[2] Pardo, L., Pardo, M. C., and Zografos, K. JJSS, 29, (1999), 213–228.
[3] Rukhin, A. L. 12th Prague Conference on Information Theory, (1994), 214–
219.
[4] Zografos, K. International Journal of Mathematics and Statistical Science 2,(1993), 5–12.

境界バイアスを回避する非対称カーネル法の様々な適用について
柿沢佳秀（北海道大学大学院経済学研究院）
１．はじめに データ {X1, . . . , Xn } ∼ iid F に基づいて，f (= F  )，F ，f (r) を推定したい．原点

∞
∞
∞
対称で 2 乗可積分なカーネル k (ただし， −∞ k(t) dt = 1， −∞ tk(t) dt = 0， −∞ t2 k(t) dt = 0を仮定する) に対して，バンド幅を h とした尺度型カーネル kh (s) = k(s/h)/h を導入する．
標準的なカーネル密度推定量 (Rosenblatt(1956)，Parzen(1962))fhRP (x) = n−1nkh (x − Xi )i=1とそれを積分または微分して得られる分布推定量 (Nadaraya(1964))FhRP (x) = x−∞fhRP (y) dy = n−1n  x−Xi
i=1−∞kh (t) dt ,密度微分推定量 (Bhattacharya(1967))
r

d
(r)RP
(x) = r fhRP (x) = n−1
fh
dx

n


(r)

kh (x − Xi )i=1の漸近特性は supp(f ) の内点 x でのみ保証される．もし supp(f ) = R なら，一般に，境界点の近くで一致性がない．実際，本報告の非負データ設定 supp(f ) = [0, ∞) において u = f (r) (orF ) を考えるとき，畳み込み積分 (kh ∗ u)(x) は ∞1 x − sku(s) ds =
h
−∞ h

 ∞
0

1 x − sku(s) ds =
h
h

 x

h

−∞

k(t)u(x − ht) dtのように変形され，非負カーネル関数 k が [−1, 1] (say) の台を持つ場合，h → 0 のとき

 min( x ,1)h−1k(t)u(x − ht) dt − u(x)

⎧
 1
⎪
1
⎪
2

⎪ h u (x)t2 k(t) dt ,x ∈ [h, ∞) ,
⎨
2
−1
≈
 1
 1
⎪
⎪

⎪
−
u(0)
k(t)
dt
+
hu
(0)
(−p + t)k(t) dt , x = hp (p ∈ [0, 1])⎩ppとなり，u(0) > 0(< 0) ならば，x = 0 の近くで深刻な O(1) の下方 (上方) バイアスがある．
一方，例外的には，u(0) = 0 のとき漸近不偏で一致性があるが，それでも，u(0) = 0，かつ，u (0) = 0 ならば, x = 0 の近くで O(h) のバイアスがある (そのオーダーは，境界点から離れた位置での O(h2 ) よりも遅い)．
この “境界バイアス問題” を対処する研究もあり，台の変換法 (Marron and Ruppert(1994))の他，境界カーネル法 (renormalization 法，reﬂection 法，一般化ジャックナイフ法は，例えば，Jones(1993) を参照) が議論された．一方，Silverman(1986; p28) は，単純な１つのアイディア
として，推定したい関数の台とマッチするような (可変的) カーネルを選ぶべきと指摘している．

有界区間 [0, 1]，半無限区間 [0, ∞) のデータに対する，このアプローチは Chen(1999,2000) がベータ/ガンマカーネル密度推定量を提案した以降，“非対称カーネル法” と知られている．

２．非負データに対する非対称カーネル法 本報告で，非負データ {X1 , . . . , Xn } ∼ iid Fに基づく f (= F  ) と F の非対称カーネル推定法について，先行研究をまずレビューし，次にf (r) の推定量の漸近特性を調べた (この構成のアイディアは先行研究と異なることに注意する)．
以下，バイアスと分散を制御する平滑化パラメータを b とする．
2.1. 非対称カーネル Chen(2000) のガンマカーネル以降，種々の非対称カーネルが個別に提案され，現在オプションは豊富になっている．本報告では，例として，対称密度 Cg g(u2 ) を土台とした非心 qBS カーネル族と qIG/qRIG カーネル族 (Kakizawa(2018,2021)) を紹介した．
2.2. 非対称カーネル密度推定量 非対称カーネル密度推定量 fb (x) = n−1 ni=1 k(Xi ; b, x)の漸近特性は Kakizawa(2021) などを参照．
2.3. 非対称カーネル分布/密度微分推定量のナイーブ版 ナイーブには，“fb (·) を積分 (or微分) すればよい” と思うに違いない，すなわち，
 x
Fb (x) =(r)0fb (y) dy = n−1

n  x

i=1

0

k(Xi ; b, y) dy ,nfb (x) =

 ∂r
dr 
−1
(x)
=
n
k(Xi ; b, x)
f
b
dxr
∂xr
i=1

という推定量もありえるが，以下の理由から，これらは推奨されない．

• 前者 (Chekkal et al.(2023)) では y 積分が数値積分にならざるをえない；• ガンマカーネルを用いた r = 1 の推定量 (Dobrovidov and Markovich(2013)) は高階微分への拡張に不利 (対数ガンマ関数の高階導関数が関与し，バイアス導出も見通しが悪い)．
 x ∞2.4. 非対称カーネル分布/密度微分推定量k(s; b, y) dy の替わりに，k(t; b, x) dt0sを使用し，Mombeni et al.(2021) は非対称カーネル分布推定量Fb (x) = n−1

n  ∞

i=1

k(t; b, x) dt

Xiを考察している．本報告は，柿沢 (2022; 連合大会/日本数学会秋季総合分科会) による非対称
カーネル密度微分推定量(r) (x) = n−1
f
b

n

(−1)r k (r) (Xi ; b, x)i=1に対して，バイアスが 2.2 節の結果で f, f  , f  を f (r) , f (r+1) , f (r+2) に読み替えて得られることを指摘した．さらに，(少なくとも) 標準正規密度を土台とした中心 qBS カーネルの場合で漸近分散も明示的に与えた．

３．おわりに 積型カーネルの使用から，多次元化は容易であるが，非積型カーネルを使うこともできる (Kakizawa(2022) を参照)．実際，qBS 型密度は，分布論として楕円密度を土台にした定式化が可能で，多変量正規密度を除く楕円密度を採用すれば，必然的に非積型となる．

Some estimators in the ADCINAR(1) process
Xiaoqiang Zeng and Yoshihide Kakizawa
Graduate School of Economics and Business Administration, Hokkaido University

1

Introduction

The analysis of count time series has rapid progress. There is a huge literature about formulationof models, probabilistic aspects, and statistical inference (see Weiß (2008,2018)). Al-Osh and Alzaid(1987) proposed nonnegative integer-valued autoregressive process of the first-order (INAR(1)), based on
the binomial thinning operator due to Steutel and van Harn (1979). Du and Li (1991) gave the stationary
condition of the 𝑝th-order INAR process, and proved its ergodicity. Silva and Oliveira (2004) discussed
some moment structures for the stationary INAR(1) process under the Poisson innovation, focusing on
Yule–Walker (YW) type equation of the third raw automoment (and autocumulant) function. Schweer
and Weiß (2016) derived, for any positve integer 𝑟, the (𝑟 + 1)th autocumulant function of the stationary
INAR(1) process under the Poisson innovation. Under a general innovation, Zeng and Kakizawa (2022)
gave the third, fourth, fifth, and sixth autocumulant functions.
Nastić et al. (2017) (see also Ristić et al. (2013)) not only introduced an alternative dependent counting
nonnegative INAR process of the first-order (ADCINAR(1)), with a modification of the original Steutel
and van Harn’s (1979) thinning operation, but also considered the estimation of the parameter 𝛼 and thenew parameter 𝜗. Recently, even for the stationary ADCINAR(1) process, Zeng and Kakizawa (2023)
studied higher autocumulant functions, except for the autocovariance function. Since the asymptotic
normality of the estimator for 𝜗 was established by Nastić et al. (2017) under an unrealistic assumption,i.e., 𝛼 is known, our first goal is to revisit such an asymptotic theory.
In the statistical analysis of the count time series data, the conditional least squares (CLS) method due
to Klimko and Nelson (1978) has been widely used. See, e.g., Al-Osh and Alzaid (1987) and Park andOh (1997) for the stationary INAR(1) process, and Nastić et al. (2017) for the stationary ADCINAR(1)process. The YW estimator for 𝛼 can be applied easily, since 𝛼 is the autocorrelation at lag 1 of various
INAR(1) type models, as in the usual stationary AR(1) process (e.g., Brockwell and Davis (1987)). Theseestimators have desirable asymptotic properties, but are biased in a finite-sample. Some authors derived
asymptotic expansions of the biases in order to perform an analytical bias-correction for the stationary
INAR(1) process (see Bourguignon and Vasconcellos (2015), Weiß and Schweer (2016), and Zeng and
Kakizawa (2022)). Our second goal is to develop the bias-corrections of the CLS and YW estimators for
𝛼 in the stationary ADCINAR(1) process.

2

ADCINAR(1) process

As usual, 𝐵( 𝑝) is the Bernoulli random variable, i.e., 𝑃[𝐵( 𝑝) = 1] = 1 − 𝑃[𝐵( 𝑝) = 0] = 𝑝(∈ [0, 1]).
Let 𝑆 𝑗 (𝛼, 𝜗) = 𝐵 𝑗 (𝜗)𝐵(𝛼/𝜗), where 0 ≤ 𝛼 ≤ 𝜗 ≤ 1(𝜗 ≠ 0) and {𝐵 𝑗 (𝜗)} is a sequence of IID Bernoulli
random variables, which is independent of 𝐵(𝛼/𝜗). An alternative generalized binomial thinning operator
was recently introduced by Nastić et al. (2017), as follows. Given a nonnegative integer-valued randomÍ
variable 𝑌 , let 𝛼⋄ 𝜗 𝑌 = 𝑌𝑗=1 𝑆 𝑗 (𝛼, 𝜗), 𝑌 = 1, 2, . . . (when 𝑌 = 0, 𝛼⋄ 𝜗 𝑌 = 0), where {𝐵 𝑗 (𝜗)} and 𝐵(𝛼/𝜗)are independent of 𝑌 . Note that 𝑆 𝑗 (𝛼, 𝜗) is distributed as the Bernoulli distribution Bin(1, 𝛼), whereas{𝑆 𝑗 (𝛼, 𝜗)} is, in general, a dependent sequence, i.e., for 𝑖 ≠ 𝑗, 𝐶𝑜𝑣 [𝑆𝑖 (𝛼, 𝜗), 𝑆 𝑗 (𝛼, 𝜗)] = 𝛼(𝜗 − 𝛼).

Nastić et al. (2017) thus defined the ADCINAR(1) process by 𝑌𝑡 = 𝛼 ⋄ 𝜗 𝑌𝑡 −1 + 𝜀 𝑡 , 𝑡 = 0, ±1, . . . , where an
innovation {𝜀 𝑡 } is a sequence of IID nonnegative integer-valued random variables, such that 𝜀 𝑡 and 𝑌𝑡 −𝑖
are independent for all integer 𝑡 and positive integer 𝑖 (it is implicitly assumed that the operations 𝛼⋄ 𝜗
at different times are performed mutually independently, which are also independent of {𝑌𝑡 } and {𝜀 𝑡 }).
Note that the INAR(1) process is a special case of the ADCINAR(1) process with 𝜗 = 𝛼(∈ (0, 1]). As
mentioned in Nastić et al. (2017), the ADCINAR(1) process {𝑌𝑡 } is strictly stationary and ergodic when0 ≤ 𝛼 ≤ 𝜗 < 1(𝜗 ≠ 0). We assumed this assumption throughout this paper.

3

Estimation of 𝛼 and 𝜗Í𝑛Í𝑛Minimizing the criterion 𝐽 (𝛼, 𝜇 𝜀 ) = 𝑡=2(𝑌𝑡 − 𝐸 [𝑌𝑡 |𝑌𝑡 −1 ]) 2 = 𝑡=2(𝑌𝑡 −𝛼𝑌𝑡 −1 − 𝜇 𝜀 ) 2 , the CLS estimatorb𝛼𝐶 𝐿𝑆 for 𝛼 can be obtained (see also Klimko and Nelson (1978) and Al-Osh and Alzaid (1987)). Further,
𝛼1,1 (see Zeng anda general estimator b
𝛼𝑐1 ,𝑐2 for 𝛼 is available, where the YW estimator is given by b
Kakizawa (2022)). The estimators b
𝛼𝐶 𝐿𝑆 and b
𝛼𝑐1 ,𝑐2 have the desirable asymptotic properties (strong
consistency and asymptotic normality). On the other hand, the criterion 𝐽 (𝛼, 𝜇 𝜀 ) is free of the parameter
𝜗, hence, another tool, referred to as the two-step CLS (2CLS) method, is needed to estimate 𝜗, as inKarlsen and Tjøstheim (1988) (see also Nastić et al. (2017)). Suppose that an estimator b
𝛼 for 𝛼(> 0)
Í𝑛
Í𝑛
−1
2
2
−1
is available. Define 𝑌 = 𝑛
𝜎𝑌 = 𝑛𝑡=1 (𝑌𝑡 − 𝑌 ) . Then, the 2CLS estimator for 𝜗
𝑡=1 𝑌𝑡 and bb
can be constructed in the closed-form, denoted by 𝜗 ( 𝛼b,𝑌 , 𝜎
b𝑌2 ) (the detail is omitted here). Under suitableb
is strongly consistent and asymptoticallyconditions, we showed that the 2CLS estimator 𝜗
2
b𝑌 )
(𝛼
b,𝑌 , 𝜎normal.

4

Bias-corrections of the CLS and YW estimators for 𝛼

Let 𝜇𝑌 (𝑢, ℓ, ℓ) = 𝐸 [(𝑌𝑡 − 𝜇𝑌 ) (𝑌𝑡+𝑢 − 𝜇𝑌 ) (𝑌𝑡+ℓ − 𝜇𝑌 ) 2 ] be fourth automoment function at lag (𝑢, ℓ, ℓ) of{𝑌𝑡 }, where 𝑢 = 0, 1 and ℓ = 1, . . . , 𝑛 − 1. Recall b
𝛼𝑌 𝑊 = b
𝛼1,1 . In the same way as Zeng and Kakizawa−2(2022), it can be shown that b
𝛼𝐶 𝐿𝑆 = b
𝛼1,0 + 𝑂 𝑝 (𝑛 ), and that the bias of the estimator b
𝛼𝑐1 ,𝑐2 is given by
Í∞
−1
4
−1
−1
𝐸 (b
𝛼𝑐1 ,𝑐2 −𝛼) = −𝑛 [1+(𝑐 1 +𝑐 2 )𝛼]−(𝑛𝜎𝑌 ) {𝑀 (1)−𝛼𝑀 (0)}+𝑜(𝑛 ), where 𝑀 (𝑢) = ℓ=1
𝜇𝑌 (𝑢, ℓ, ℓ),𝑢 = 0, 1. We proceeded to estimate 𝑀 (1) − 𝛼𝑀 (0) in two different ways. One is based on the analyticalevaluation; the other is the nonparametric method applying the lag window-type kernel (note that theformer analytical evaluation needs a further extra tedious calculation of 𝑀 (1) − 𝛼𝑀 (0)).

5Conclusion and future worksb
The asymptotic normality of the 2CLS estimator 𝜗b𝑌2 ) for 𝜗 was revisited. We proposed the lag(𝛼b,𝑌 , 𝜎window-type bias-correction and the analytical bias-correction of the commonly used CLS and YW
estimators for 𝛼 in the stationary ADCINAR(1) process and then conducted some simulations.
Bootstrap and jackknife procedures are alternative methods of constructing bias-corrected estimators.
Jentsch and Weiß (2019) proved the validity of bootstrapping for the stationary INAR process, so that theextension of thier result to the stationary ADCINAR(1) process was left for a future topic.

Quantification and Formulation of Galaxy Evolution by
Manifold Learning
Tsutomu T. Takeuchi1,2∗, Suchetha Cooray3†, Hai-Xia Ma1 , Kiyoaki Christopher Omori1‡,
Wen Shi1 , Ryusei R. Kano1 , Sena A. Matsui1 , and Aina May So4,1
1. Division of Particle and Astrophysical Science, Nagoya University, Japan,
2. Research Center for Statistical Machine Learning, Institute of Statistical Mathematics, Japan,
3. National Astronomical Observatory of Japan,
4. Department of Physics, Gakushuin University, Japan

Matter in the early Universe was almost uniform. Slightly dense regions have grown by gravity, and finallyturned into galaxies. Astrophysicists have searched for the equations governing the complex physical phenomenaof galaxy formation and evolution over 13 billion years of the cosmic age, from the first principles of physics.
We, instead, elucidate the physics of galaxy evolution by applying manifold learning, one of the latest methodsof data science, to a feature space spanned by galaxy luminosities and cosmic time. We have discovered thatgalaxies consist a complex but low-dimensional subset in an ultra-high-dimensional feature space (Siudek et al.,2018). We refer to it as the galaxy manifold in a modern sense. The discovered galaxy manifold by Siudek
et al. (2018) has a strongly nonlinear one-dimensional structure. Further, surprisingly, the spectra of the sample
galaxies were discriminated only by a few broadband luminosities, not complicated combinations of quantities.
This suggests that the galaxy evolution at optical wavelengths can be described by only a few parameters at
most. This is a new characterization of galaxy evolution that could have never been found by conventional
methods. We will further pursue the galaxy manifold and elucidate the dependence of parameters (probably
several at most) that govern the physics of galaxy evolution from the galaxy manifold, and derive the governing
equations for galaxy evolution.
We then applied the unsupervised machine learning to express it even in an easier way. This is a part of
the technique referred to as the dimensionality reduction. We used the method topological data analysis toelucidate the (possibly several at most) parameter dependence that governs the physics of galaxy evolution
from galaxy manifolds. Especially, a method known as the manifold learning is optimal for parametrizing thegalaxy manifold. Multimodality, and dispersions in classical scaling laws are merely a consequence of unsuitableprojection of the galaxy manifold. We make use of this method to construct the governing equation(s) of galaxy
evolution.
After obtaining the galaxy manifold, for example, observational quantities as luminosities missing for objects,
or more physical quantities as SFR and stellar mass M∗ can be straightforwardly estimated from the locationon the manifold. We adopt the algorithm Isomap Tenenbaum et al. (2000) and UMAP (Uniform Manifold
Approximation and Projection) McInnes et al. (2018, 2020). Isomap defines the neighboring points by using
input-space distance, and the distant points as a sequence of “short hops” between neighboring points. Isomap
tries to find shortest paths in a graph with edges connecting neighboring data points By construction, Isomap
preserves the “surface density” of data points in the feature space. UMAP is based on differential geometry
and algebraic topology. The algorithm is founded on three assumptions: 1) the data are uniformly distributedon a Riemannian manifold, 2) the Riemannian metric is locally constant (or can be approximated as such),∗ E-mail: tsutomu.takeuchi.ttt@gmail.com.
† JSPS Fellow (PD)‡ JSPS Fellow (DC2)1and 3) the manifold is locally connected. From these assumptions it is possible to model the manifold with afuzzy topological structure. Sine it defines the manifold so that the data points distribute as homogeneously as
possible, it does not preserve the surface density of data points. UMAP also preserves some important structuralproperties, and it is more robust against noise than Isomap. Manifold learning algorithm can “unfold” a curved
and/or rolled manifold in the feature space, and provide a local coordinate system on it Tenenbaum et al. (2000);
Roweis & Saul (2000). We also stress that two different algorithms, Isomap and UMAP yield similar manifolds.
Since Isomap preserves the density of data point cloud, we observe that the manifold has a density structure,i.e., dense and sparse regions on the manifold.
To demonstrate that the manifolds contain important information on galaxy evolution, we compared SFR
and stellar mass as a function on the manifold. The behaviors of SFR and stellar mass are both qualitatively
very similar, suggesting that the estimated manifold structure is robust. This means that the manifold learning
actually “learns” the important characteristics of galaxy evolution at optical wavelengths. Since we could connect
the manifold and some physical quantities such as SFR and stellar mass, we can also parametrize the evolutionof galaxies on the manifold. As stars are formed, the stellar mass, namely the total mass of accumulated stars
will increase. This is one of the fundamental aspect of galaxy evolution, and we can visualize this evolution asa vector field on the manifold.
Last question is how to describe and interpret the evolutionary track of galaxies on the manifold. We applieda classical theoretical model of the chemical evolution of galaxies. Chemical evolution is a field of galactic physics
that treats the formation and evolution of elements in galaxies based on the stellar evolution theory. This strange
terminology comes from the fact that the theory has been used to analyze the chemical composition of stars andthe interstellar medium (ISM) in a galaxy. The important physical process is the nucleosynthesis via nuclear
fusion in the center of stars. We adopt a simple model with infall and outflow of matter proposed by Lilly et al.
(2013).
M∗ (tn+1 ) = M∗ (tn ) + (1 − r)SFR(tn )∆t ,(1)

Mgas (tn+1 ) = Mgas (tn ) − (1 − r + η)SFR(tn )∆t(2)where r is the returning mass fraction, η is the mass loading factor, and ∆t is the timestep (Cooray et al., 2023).
We calculate theoretical evolutionary tracks of galaxies from eqs. (1) and (2). Interpretation of the vector fieldon the galaxy manifold can be given by comparing eqs. (1) and (2) with the vector field on the galaxy manifold.
However, a sophisticated methodology is desired to interpret the observed manifold more directly.
We restricted our discussion to multiwavelength photometric surveys, but this method can in principle beextended to spectroscopic surveys. In this work, we have discussed the optical properties, but also we can include
any other properties of galaxies, not only luminosities at other wavelengths but also dynamical, structural, andenvironmental properties. Manifold learning will provide a fundamentally new insight to the studies on theformation and evolution of galaxies.

References
Cooray, S., Takeuchi, T. T., Kashino, D., et al. 2023, Mon. Not. R. Astron. Soc. , 524, 4976
Lilly, S. J., Carollo, C. M., Pipino, A., Renzini, A., & Peng, Y. 2013, Astrophys. J. , 772, 119
McInnes, L., Healy, J., & Melville, J. 2020, UMAP: Uniform Manifold Approximation and Projection for
Dimension Reduction, arXiv:1802.03426
McInnes, L., Healy, J., Saul, N., & Großberger, L. 2018, Journal of Open Source Software, 3, 861
Roweis, S. T., & Saul, L. K. 2000, Science, 290, 2323
Siudek, M., Malek, K., Pollo, A., et al. 2018, Astronomy and Astrophysics, 617, A70
Tenenbaum, J. B., Silva, V. d., & Langford, J. C. 2000, Science, 290, 2319
2

GMANOVA モデルにおけるフルランク仮定が不要な罰則無推定法(報告書)

中京大学 教養教育研究院 永井 勇
本講演では n 個の各個体に対して，時間と共に p 回測定されたデータは経時測定データの分析を考えていた．このような経時測定データにおいて，全ての個体で測定した時点が揃っている場合，Potthoﬀ and Roy (1964) により提案された次の一般化多変量分散分析（GMANOVA）モデルでの分析がよくされていることを紹介した;

Y = 1n µ′ X ′ + AΞX ′ + E ，(1)

ここで，Y は各行が各個体の経時測定データからなる n × p 既知行列，1n は全ての要素が

1 の n 次元のベクトル，µ は q 次元未知ベクトル，X は個体内計画行列と呼ばれる測定時点から作られる p × q 既知行列 (詳細は後述)，A は個体間計画行列と呼ばれる測定時点に無関係な各個体の特徴を表す変数からなる n × k 既知行列で各列で中心化されているとし，

Ξ は k × q 未知行列であり，E は E[E] = 0n 0′p で Var[vec(E)] = Σ ⊗ In の n × p 誤差行列であり，0r は全ての要素が 0 の r 次元ベクトル，Σ は未知の p × p 正則行列である．

Y の経時的な変動 (経時変動) を推定する際には，個体内計画行列 X の i 行目として，i番目の測定時点 ti (i = 1, . . . , p, t1 < · · · < tp ) の関数からなる行列が用いられることを紹
q−1

介した．特に，本講演では分かりやすさのために，X の i 行目として (t0i , . . . , ti

) を用いて未知の µ や Ξ を推定することで，切片と各次数の係数に対応した部分が推定でき，経時変動を ti の (q − 1) 次多項式で推定できることを例示した．
このモデル (1) において，µ や Ξ の推定は次のリスクを最小にすることで行われることを紹介した;

tr{(Y − 1n µ′ X ′ − AΞX ′ )Σ−1 (Y − 1n µ′ X ′ − AΞX ′ )′ }．

(2)

このリスクを最小にする µ や Ξ を求める際に，従来は rank(A) = k かつ rank(X) = qを仮定していることを紹介し，本講演ではこれらの仮定が満たされない場合の推定法を提案した．そこでは，永井 (2021) のアイデアを用いて，シンプルに A や X を (A1 , . . . , Ad1 )や (X1 , . . . , Xd2 ) と分割し，µ や Ξ をそれぞれの大きさに合わせて (µ1 , . . . , µd2 ) や Ξ の j
行目を (Ξj1 , . . . , Ξjd2 ) (j = 1, . . . , d1 ) と分割することを考えた．この分割した行列などを用いて，µ や Ξ を推定する代わりに µ1 , . . . , µd2 や Ξ11 , . . . , Ξ1d2 , . . . , Ξd1 1 , . . . , Ξd1 d2 を推定することで，rank(A) = k かつ rank(X) = q より非常に緩い仮定の下で，罰則を付けずにリスク (2) を最小にする µ や Ξ の推定法を提案した．さらに，分割した行列を用いてリスクを直接最小にする Ξij (i = 1, . . . , d1 ; j = 1, . . . , d2 ) などを求めると，それらを得るた
めに求めようとしている部分以外の全ての部分が必要となり，初期値や反復計算が必要となることを紹介した．この問題を解決するため，リスクの展開を工夫した形を構築し，それに基づいた推定法を提案した．

引用文献:[1] Katayama, S. & Imori, S. (2014) Lasso penalized model selection criteria for highdimensional multivariate linear regression analysis. J. Multivariate Anal., 132, 138–
150.
[2] 永井 勇 (2021) 高次元小標本における多変量線形回帰モデルでの推定法，2021 年度統計関連学会連合大会 講演.
[3] Potthoﬀ, R. F. & Roy, S. N. (1964) A generalized multivariate analysis of variancemodel usefule specially for growth curve problems. Biometrika, 51, 313–326.
[4] Srivastava, M. S. (2002) Methods of Multivariate Statistics, Wiley-Interscience.
[5] Yanagihara, H., Nagai, I., Fukui, K. & Hijikawa, Y. (2023) Modiﬁed Cp criterion
in widely applicable models. Smart Innov. Syst. Tec., 352, Intelligent Decision
Technologies 2023: Proceedings of the 15th KES International Conference on Intelligent Decision Technologies KES-IDT-2023 (eds. I. Czarnowski, R. J. Howlett & L.
C. Jain), 173-182.

Two step estimations via the Dantzig selector for models of
stochastic processes with high-dimensional parameters
Kou Fujimori1 and Koji Tsukuda2
1

Faculty of Economics and Law, Shinshu University.
2
Faculty of Mathematics, Kyushu University.

Let θ ∈ Θ ⊂ Rp be an unknown parameter of interest and h ∈ H a possibly inﬁnitedimensional nuisance parameter where H is a metric space equipped with a metric dH .
Consider the following random maps:
ψn(1) : Θ → Rp ,ψ̃n(1) : Θ → Rp ,

Ψ n : Θ × H → Rp ,Ψ̃n : Θ × H → Rp ,(1)

where n is a number of observations. The random maps ψn and Ψn are corresponding to(1)some score functions and ψ̃n and Ψ̃n are their compensators, respectively. We supposethat the true values θ0 ∈ Θ and h0 ∈ H satisfy thatψ̃n(1) (θ0 ) ≈ 0,

Ψ̃n (θ0 , h0 ) ≈ 0.

We are interested in the estimation problem for θ0 = (θ01 , . . . , θ0p )⊤ under the following
high-dimensional and sparse settings;
We ﬁrst construct an estimator for T0 to achieve the dimension reduction. To do this,(1)(1)
we consider the following Dantzig selector type estimator θ̂n for θ0 , based on ψn :θ̂n(1) := arg min ∥θ∥1 ,θ∈Cn

Cn = {θ ∈ Θ : ∥ψn(1) (θ)∥∞ ≤ λn },(1)where λn is a tuning parameter. Using θ̂n , we deﬁne T̂n as follows:T̂n := {j : |θ̂nj | > τn },where τn is a threshold level, which is a tuning parameter. Then, we have the following
theorem.

1

Theorem 1. Under some regularity conditions, it holds that
(
)
2
(1)
P ∥θ̂n − θ0 ∥∞ > cn → 0, n → ∞,
δ
(1)

(1)

(0.1)

(1)

where cn := ∥ψn (θ̂n ) − ψn (θ0 )∥∞ . Especially, it holds that∥θ̂n(1) − θ0 ∥∞ = Op (λn ),n → ∞.

Moreover, if the threshold τn satisﬁes that 4λn /δ < τn < inf j∈T0 |θ0j |/2, then it holds that
P (T̂n = T0 ) → 1,n → ∞.

(0.2)

For every index set T , we consider the following random map restricted by T :ΨnT : ΘT × H → R|T | ,Ψ̃nT : ΘT × H → R|T | ,where ΘT is a set of sub-vectors of Θ restricted by T . Let ĥn be an estimator of h ∈ Hsuch that dH (ĥn , h0 ) = op (1) as n → ∞. Then, we consider the new estimator θ̃n for theparameter θ of interest, with help of ĥn and T̂n , as a solution to the following equations:ΨnT̂n (θ̃nT̂n , ĥn ) ≈ 0,

θ̃nT̂ c = 0.
nThen, we can derive the asymptotic distribution of θ̃n , as well as Theorem 2.1 of Nishiyama(2009). In this talk, we apply the general theory discussed above to ergodic time seriesmodels and models of diﬀusion processes.

Acknowledgements
This work was supported by JSPS KAKENHI Grant Numbers 21K13271(K.F.), 21K13836(K.T.).

ReferencesBühlmann, P. and van de Geer, S. A. (2011). Statistics for high-dimensional data: methods, theory and applications. Springer Science & Business Media.
Candes, E. and Tao, T. (2007). The Dantzig selector: statistical estimation when p ismuch larger than n. Ann. Statist. 35 2313–2351.
Huang, J., Sun, T., Ying, Z., Yu, Y., and Zhang, C.-H. (2013). Oracle inequalities for the
LASSO in the Cox model. Ann. Statist. 41 1142–1165.
Nishiyama, Y. (2009). Asymptotic theory of semiparametric Z-estimators for stochastic
processes with applications to ergodic diﬀusions and time series. Ann. Statist. 37
3555–3579.
Tibshirani, R. (1996). Regression shrinkage and selection via the lasso. J. Roy. Statist.
Soc. Ser. B 58 267–288.
2

Statistical mathematics of biodiversity indices for
spatial species distribution model
Shimatani, Ichiro Ken (The Institute of Statistical Mathematics)
In order to express biodiversity, many indices were proposed. The number of different
species (species richness) has been most widely used, and when we have species
abundance data, the following two indices have been most commonly used.
Shannon entropy H = −Σs ps log ps ,Simson index D = 1 − Σs ps2 ,where ps is the relative frequency of species s (=the probability that a randomly chosen
individual is species s).
It is shown that the three indices are unitedly interpreted in the Hill number defined by1Sq

D = ( ps ) 1−qqs =1When q = 0, it is identical to the species richness. When q = 1, it is equivalent to the
Shannon entropy;1−  p s ln p s
D=e s,

and when q = 2, it corresponds to the Simpson index;
1
D=
.
p
2

S

2

s =1

s

Thus, we can comparatively examine biodiversity of two or more communities not by asingle number but by curves parametrized by 0 ≤ q (practically, in many case, ≤ 2 or 3),
which indicates a balance (weight) between rare species and dominant ones.
The Hill number, however, does not provide a perfect quantification of biodiversity. The
biggest problem is that the Hill number is highly underestimated if observed relative
frequencies are inserted into the equation (Fig. 1).
Other indices are also proposed and comparatively examined. An example is thegeneralized Simpson’s entropy defined as;S r =  ps (1 − ps ) r(r = 1, 2, …)

s =1

This is interpreted as the probability that (r + 1)st observation will be of a species that has
not been observed before, and an unbiased estimator is known.

Fig. 1 Hill number is highly biased if observed relative frequencies are directly used in the equation. Examples using
100 samples. The left panels indicate the given relative frequencies and the Hill numbers for the four ecologicalcommunities. The upper right two graphs are for community A (left) and B (right). The lower two are for community
C (left) and D (right).
50

Estimates of Hill number

1.0

Relative frequency

0.8

0.6

50
A TrueA meanA 1.96SD
A 1.96SD

40

30

40

30

20

20

10

10

0.4

0

0
0

0.2

0.0
B
C
Community

140

A

B

C

D

120
100
Hill number

1

1.5

2

50

50

40

40

30

30

20

20

10

10

0

0.5

1

1.5

2

0

0.5

1
q

1.5

2

D

Estimates of Hill number

A

0.5

80
60
40

0

20

0
0

0.5

0
0

0.5

q 1

1.5

1
q

1.5

2

2

Recently. a goal of community ecology tends to change from finding the “true”
community to identifying a hidden stochastic process that controls realized communities
and an observed community is nothing more than one realization of the stochastic process(Fig. 2).

Stochastic process

Fig. 2e.g. spatial point process
Realization of point patterns

Given byequation

Diversity
index
Classical
diversity studies(e.g. unbiased
estimator)

Real community (true abundances ofmember species, N1, …, NS)

Observation model(sampling)

Classical
statisticalinference

Up-to-datestatisticalinference

Observed abundances(numbers of individuals, n1, …, nS)

In this regard, the presence of an unbiased estimator may not be so big advantage and a
new overall statistical inference is an on-going statistical issue in biodiversity studies..

Simultaneous estimation of Poisson means: recent developments and its applications
Yuan-Tsung Chang (Mejiro University), Shinozaki Nobuo (Keio University)
In estimating p(≥ 2) independent Poisson means, Clevenson and Zidek (1975) have proposed a class of estimatorsthat shrink the unbiased estimator to the origin and dominate the unbiased one under the normalized squared
error loss. This class of estimators was subsequently enlarged in several directions. Here, we disscuss
1. the problem and proposes new classes of dominating estimators using prior information pertinently. Dominance is shown by partitioning the sample space into disjoint subsets and averaging the loss difference overeach subset. Estimation of several Poisson mean vectors is also discussed. Further, simultaneous estimationof Poisson means under order restriction is treated and estimators which dominate the isotonic regression
estimator are proposed for some types of order restrictions.(Chang and Shinozaki (2019)).
2. the shrinkage estimation of Poisson means when observations are given in the form of a two-way contingency
table. Assuming a multiplicative Poisson model, estimators which shrink to the specified values or an order
statistic in one dimension and in two dimensions are considered and are shown to dominate the maximum
likelihood estimator (MLE) under normalized squared error loss. Further, assuming the full model, shrinkage
to the multiplicative model is devised to improve upon the unbiased estimator. Shrinkage is made after
determining the basic cells so that the observed frequency is not smaller than the estimated frequency foreach of the other cells. (Chang and Shinozaki (2022)).

1

Some classes of improved estimators using prior information

Let Xi be distributed as P o(λi ), i = 1, . . . , p, and suppose that X1 , . . . , Xp are statistically independent. Let theloss function when we estimate λ = (λ1 , . . . , λp ) by λ̂ = (λ̂1 , . . . , λ̂p ) be the normalized squared error one
L(λ, λ̂) =pX(λ̂i − λi )2
i=1

λi

.

(1.1)

Then Clevenson and Zidek (1975) were the first to propose a class of estimators of the formλ̂CZi (X) = Xi −where Z =Pp

i=1 Xi .

ϕ(Z)Xi ,Z +p−1i = 1, . . . , p,

(1.2)

Clevenson and Zidek (1975) have shown that when p ≥ 2 and if ϕ(·) is a non-decreasing
CZt
CZ
function satisfying 0 ≤ ϕ(·) ≤ 2(p − 1), then λ̂ (X) = (λ̂CZ
1 (X), . . . , λ̂p (X)) dominates X. Since then broadclasses of dominating estimators have been given by many authors, including Tsui and Press (1982), Hwang (1982),
Ghosh et al. (1983), and Chou (1991).

1.1

Shrinking to a specified point other than the origin

We first consider a class of estimators which shrink Xi to ai only when xi ≥ ai , i = 1, . . . , p, where ai 0 s are specified
nonnegative values which are chosen according to prior information about λi 0 s. Let C = {(x1 , . . . , xp ) | xi ≥ ai , i =1, . . . , p} and let IC be its indicator function. We consider a class of estimators of the form
λ̂ai (X) = Xi − ϕ(Z)Xi − aiIC ,Z +di = 1, . . . , p,

Pp
where Z = i=1 (Xi − ai ) and d is a positive constant.
Theorem 1.1. Let ϕ(·) be a non-decreasing function which satisfies 0 ≤ ϕ(·) ≤ 2(p − 1) and suppose thatad ≥ sup ϕ(·)/2. Then λ̂ (X), whose i-th component is given by (2.1), dominates X under the normalized squared
error loss (1.1).

1.2

Shrinking to the order statistics

Here we consider a class of estimators which shrink to the order statistics. We first consider the following shrinkage
estimator toward X(1) = min{X1 , . . . , Xp } :(1)λ̂i (X) = Xi − ϕ(W )Xi − X(1),W +d1

i = 1, 2, . . . , p,Ppwhere W = i=1 (Xi − X(1) ). Then we have the following.
Theorem 1.2. Suppose that p ≥ 3. If ϕ(·) is a non-decreasing function which satisfies 0 ≤ ϕ(·) ≤ 2(p − 2) andd ≥ sup ϕ(·)/2, then λ̂error loss (1.1).

2(1)(X), whose i-th component is given by (2.4), dominates X under the normalized squared

Shrinkage estimators of multiplicative Poisson means in two-way
contingency tables

We consider two-way multiplicative model where xij , i = 1, . . . , I, j = 1, . . . , J, are independent random Poisson
random variables with means
λij = λαi βj , i = 1, . . . , I, j = 1, . . . , J,PI
PJ
where αi ≥ 0 and βj ≥ 0 satisfy i=1 αi = 1 and j=1 βj = 1, respectively. We denote the one-dimensionalfrequencies and the total frequency by
xi+ =

J
X

xij , i = 1, . . . , I, x+j =

j=1

I
X

xij , j = 1, . . . , J, x++ =

i=1

I X
J
X

xij .

i=1 j=1

As discussed in Hara and Takemura (2006) complete sufficient statistics are x1 = (x1+ , . . . , xI+ ) and x2 =(x+1 , . . . , x+J ). The MLE of λij is
L
λ̂M
=
ij

( xi+ x+j
x++
0

if

x++ 6= 0

if

x++ = 0.

They have given a class of improved estimators which shrink the MLE toward the origin under the normalized
squared error loss. The simple one is)(xi+ x+jdHT, i = 1, . . . , I, j = 1, . . . , J,
δij =
1−
x++
x++ + d

2.1

One-dimensional shrinkage to an order statistic or a specified point

Let x(`)+ be the `-th smallest observation among x1+ , . . . , xI+ . We assume that I ≥ `+2 and consider the following
estimator which shrinks xi+ toward x(`)+ when xi+ ≥ x(`)+ :()(xi+ − x(`)+ )+
x+j
(1)
δij =
xi+ − ϕ(W ), i = 1, . . . , I, j = 1, . . . , J,
x++
W +d
PI
where W = i=1 (xi+ − x(`)+ )+ , a+ = max(0, a) and d is a positive constant. Then we have the following.
Theorem 2.1. Suppose that ϕ(W ) is a non-decreasing function satisfying 0 ≤ ϕ(W ) ≤ 2(I − ` − 1) and that d ≥PI(1)L
sup ϕ(W )/2. Then δij , i = 1, . . . , I improves upon the MLE λM
ij , i = 1, . . . , I under the loss function
i=1 (λ̂ij −λij )2 /λij for any j = 1, . . . , J.

2.2

One-dimensional shrinkage to a specified point

Let bi ≥ 0, i = 1, . . . , I be given numbers and we propose the following shrinkage estimator which shrinks xi+ to bi
when xi+ ≥ bi :()(xi+ − bi )+x+j(2)
xi+ − ϕ(N, W ), i = 1, . . . , I, j = 1, . . . , J,
δij =
x++
W + d(N )PI
where W = i=1 (xi+ − bi )+ and N = #{i|xi+ ≥ bi }. Then we have the following.
Theorem 2.2, Suppose that ϕ(N, W ) is a non-decreasing function of W and satisfies 0 ≤ ϕ(N, W ) ≤ 2(N − 1)+(2)
for any 0 ≤ N ≤ I. Suppose that d(N ) ≥ supW ϕ(N, W )/2. Then δij , i = 1, . . . , I improves upon the MLE
P
I
L
2
λ̂M
ij , i = 1, . . . , I under the loss function
i=1 (λ̂ij − λij ) /λij for any j = 1, . . . , J.

2

Fredholm 行列式に関する Darling の公式とその拡張
田中 勝人（一橋大学名誉教授）

Brown 運動の 2 次汎関数で表現される統計量の特性関数を導出するにあたり，特定の核関数に対しては，従来の複雑な手続きを回避して簡単で効率的な方法が可能であることを証明し，具体例を示した．内容は以下の通りである．
ほぼ定符号の核関数 K(s, t) と Brown 運動 {W (t)} に対して，統計量∫ 1∫ 1S=

K(s, t) dW (s) dW (t)00を考える．このとき，次の定理が成立する．
定理 1：Anderson-Darling の定理. 統計量 S の特性関数は，(

E eiθS

)

[

{

=

∞
∏

(

n=1

1−

}]

∫ 1∫ 1K(s, t) dW (s) dW (t)= E exp iθ
0

2iθ
λn

0

)−1/2

(

)−1/2

= D(2iθ)で与えられる．ここで，λ1 , λ2 , , . . . , は K(s, t) の固有値で重複度の分だけ繰り返される．
また，D(λ) は K(s, t) の Fredholm 行列式 (FD) である．

FD の定義から FD を導出することは一般に困難である．そのために積分方程式を同値な微分方程式（＋境界条件）に変換して，その解から FD の候補を見つけ，実際に FDであることを証明するのが通常の方法である．しかし，核関数が複雑な場合は，導出が困難となる．
ある条件をみたす核関数に限定した上で，微分方程式の議論などを回避することにより，簡単に効率的に FD を計算できる方法を提案した．そのための方法の一つとしてDarling (1955) の公式を紹介した．
定理 2：Darling の公式． 核関数

K1 (s, t) = min(s, t) − st − ψ(s)ψ(t),ψ(0) = 0,

ψ(1) = 0を定義する．ここで，ψ(t) は連続微分可能である．このとき，K1 (s, t) の FD は次式で与
えられる．
√
)
∫ 1 (∫ t√√sin λ
D1 (λ) = √−2
ψ ′ (s) cos λs ds ψ ′ (t) cos λ(1 − t) dt.
(1)
0
0
λ

Darling の公式の変形バージョンとして，次の公式を証明した．
定理 3. 核関数

K2 (s, t) = 1 − max(s, t) − ψ(s)ψ(t),ψ(1) = 0

を定義する．ここで，ψ(t) は連続微分可能である．このとき，K2 (s, t) の FD は次式で与えられる．

D2 (λ) = cos)√√√ ∫ 1 (∫ t ′√ψ (s) sin λs ds ψ ′ (t) cos λ(1 − t) dt.
λ+2 λ
0

(2)

0

さらに，核関数が次の場合の FD に関する定理を証明した．
定理 4. 核関数

K3 (s, t) = min(s, t) − st −m∑ψj (s)ψj (t),ψj (0) = ψj (1) = 0 (j = 1, . . . , m)

j=1の FD は√sin λ
D3 (λ) = √|P (λ)|(3)λで与えられる．ここで，P (λ) は m × m の対称行列で，その (j, k) 要素は次のように定
義される．
√
(∫ t
)
√
√
2 λ ∫1′√ψj (s) cos λs ds ψj′ (t) cos λ(1 − t) dt,
Pjj (λ) = 1 −0
sin λ 0√∫ 1∫ 1λPjk (λ) = − √ψj′ (s)ψk′ (t)L1 (s, t) ds dt(j ̸= k),sin λ 0 0√√
 cos λs cos λ(1 − t) (s ≤ t)
L1 (s, t) = √√cos λt cos λ(1 − s) (s ≥ t).
核関数が複雑となったもう一つの場合についても証明した．
定理 5. 核関数

K4 (s, t) = 1 − max(s, t) −m∑ψj (s)ψj (t),ψj (1) = 0

j=1の FD はD4 (λ) = cos√
λ |P (λ)|

で与えられる．ここで，√ ∫ (∫
)
1
t
√
√
2 λ
′
√
Pjj (λ) = 1 +ψj (s) sin λs ds ψj′ (t) cos λ(1 − t) dt,
0cos λ 0√∫ 1∫ 1λ√(j ̸= k),Pjk (λ) =ψj′ (s)ψk′ (t)L2 (s, t) ds dt
0
0
cos λ
√
√


 sin λs cos λ(1 − t) (s ≤ t)
L2 (s, t) = √√sin λt cos λ(1 − s) (s ≥ t).

(4)

Comprehensive Interval-Valued Time Series Model with Application to the
S&P 500 Index and PM2.5 Level Data Analysis
Liang-Ching Lina *, Hao Sunga and Sangyeol Leebab

1

Department of Statistics, Institute of Data Science, National Cheng Kung University, Tainan, Taiwan.

Department of Statistics, Seoul National University, Seoul, 08826, South Korea.

Introduction

In recent years, modern data architecture has been recorded in short time intervals, but stillcontains rich information resources. Such data may contain several problems such as missing
data, repeated observations, and recording non-equidistant time points. To resolve such issues,
researchers have considered reorganizing the data into interval-valued data consisting of daily
or weekly maximum and minimum values as it can preserve more information than merely
aggregated daily or weekly mean data. For instance, the PM2.5 levels are recorded hourly but
there are 5588 missing observations. To handle missing data, Gao and Tsay (2019) aggregatedthe data into weekly mean observations. However, we arranged the data into weekly maximum
and minimum observations. Then, we consider to fit the auto-interval-regressive (AIR) model
proposed by Lin et al. (2021) and the auto-interval-regressive moving average (AIRMA) model
by Lin et al. (2023). Further, we fit the GHVAIRMA (generalized heteroscedastic volatility
AIRMA) model to improve the accuracy of the one-step-ahead prediction of the high andlow prices of the S&P 500 index. Overall, our findings strongly confirm the adequacy of the
proposed model.

2Interval-Valued Time Series Model

Let Yt = (Yu,t , Yl,t )T . Lin et al. (2021) proposed the AIR(p) model as follows:
Yt = ϕ1 Yt−1 + · · · + ϕp Yt−p + At ,where ϕi , i = 1, . . . , p, are model parameters, At = (Au,t , Al,t )T are i.i.d. errors with Au,t =max{A1,t , . . . , An,t } and Al,t = min{A1,t , . . . , An,t }; and Ai,t , i = 1, . . . , n, are i.i.d. normal
random variables with a mean of zero and variance of σ 2 for some n > 1. Lin et al. (2023)
proposed the IVMA(q) if the following equation is satisfied:
Yt = At − θ1 At−1 − · · · − θq At−q .

(2.1){Yt } is said to follow an AIRMA(p, q) if the following equation holds:
Yt = ϕ1 Yt−1 + · · · + ϕp Yt−p + At − θ1 At−1 − · · · − θq At−q .

(2.2)

If the error terms Ai,t , i = 1, . . . , n have a mean of zero and the conditional heteroscedastic
variance of σt2 . Then, At follows a general heteroscedastic volatility model with lags r and s,* E-mail address: lclin@ncku.edu.tw1abbreviated as GHV(r, s), if σt2 satisfies equation
σt2 = α0 +rX2
αi γt−i
+

i=1

s
X

2
βj σt−jwhere γt2 =

j=1

Φ−1

Au,t − Al,t

n−α
− Φ−1
n−2α+1!21−αn−2α+1



.

(2.3)

Combined with the AIRMA model from (2.2), we have the GHV(r, s)-AIRMA(p, q), where theinterval-valued time series {Yt } satisfies
Yt = ϕ1 Yt−1 + · · · + ϕp Yt−p + At − θ1 At−1 − · · · − θq At−q .

3

(2.4)

Real Data Analysis

In this application, we consider the prediction of one-step high and low prices for the S&P 500
index. The MDEs of the proposed GHV(1,1)-IVMA(3) model outperform the VAR-DCC andHVAIR(1,1) models.
For the PM2.5 level data, we fit the AIR(p1 ), IVMA(q1 ), and AIRMA(p2 , q2 ) models withp1 ≤ 10, q1 ≤ 12, p2 ≤ 2, q2 ≤ 4, and the exogenous parameter n = 24 × 7 = 168 to the dataset,and choose the best models as those with the lowest AIC and BIC values. We investigate
the relationships among the stations by using the dandelion plot proposed by Zhang and Lin(2022), which establishes the pairwise correlations of the center-to-center and range-to-rangeof the residuals. The residuals deduced from the proposed models provide more intuitively
appealing results. For comparison, we also fit a VAR(1) model to the weekly mean data. Allother stations exhibit higher correlations with each other. The results indicate that too muchinformation has been aggregated into a single point. These comparative results indicate that
our method preserves more information and provides a more intuitively appealing result.

4

Concluding Remarks

Our models can be relaxed by replacing it with other distributions, and the coefficients of the
maximum and minimum values can be allowed to differ from each other as follows:
  
 
 

Yu,t
ϕ1 Yu,t−1Au,tθ1 Au,t
=
+
−
.
Yl,t
ϕ2 Yl,t−1Al,tθ2 Al,t−1
All these issues are worth further investigation and should be examined in future studies.

Reference
1. Gao, Z., and Tsay, R. S. (2019). A Structural-Factor Approach to Modeling HighDimensional Time Series and Space-Time Data. Journal of Time Series Analysis, 40(3),343-362.
2. Lin, L.-C., Chien, H.-L., and Lee, S. (2021). Symbolic interval-valued data analysis fortime series based on auto-interval-regressive models. Statistical Methods & Applications,30, 295-315.
3. Lin, L.-C., Hao, S., and Lee, S. (2023). Comprehensive interval-valued time series modelwith application to the S&P 500 index and PM2.5 level data analysis. Applied Stochastic
Models in Business and Industry, 39(2), 198-218.
4. Zhang, M., and Lin, D. K. (2022). Visualization for Interval Data. Journal of Computational and Graphical Statistics, DOI: 10.1080/10618600.2022.2066678.
2

Rank-R Matrix Autoregressive Models for Modeling
Spatio-Temporal Data
Nan-Jung Hsu, Institute of Statistics, National Tsing-Hua University
Hsin-Cheng Huang, Institute of Statistical Science, Academia Sinica
Ruey S. Tsay, Booth School of Business, University of Chicago
Tzu-Chieh Kao, Institute of Statistics, National Tsing-Hua University
September 30, 2023
Abstract
We develop a matrix-variate autoregressive (MAR) model to analyze spatio-temporal data organized on a regular grid in space. The model is an extension of the bilinear MAR spatial model of
Hsu, Huang and Tsay (2021) by increasing its flexibility and applicability in empirical applications.
Specifically, we propose to model each autoregressive (AR) coefficient matrix of the MAR model by R
bilinear terms, thereby establishing a rank-R model. The extension can be interpreted as decomposing the AR dynamics of the data into R bilinear MAR components. We further incorporate a bandedneighborhood structure for AR coefficient matrices and utilize a flexible nonstationary low-rank covariance model for the spatial innovation process, leading to a parsimonious model without sacrificing
its flexibility. We estimate all parameters of the model by the maximum likelihood method and develop a computationally efficient alternating direction method of multipliers algorithm, involving
only closed-form expressions in all steps. Applications to a wind-speed dataset and an employment
dataset, as well as two simulation experiments, demonstrate the effectiveness of the proposed methodin estimation, model selection, and prediction.

Keywords: Alternating direction method of multipliers, Bayesian information criterion, Kronecker product, Low-rank approximation, Matrix-variate time series, Nonstationary spatial model, Singular value
decomposition.

1

Matrix Autoregressive Spatio-Temporal Models

Consider a zero-mean n × m matrix-variate time series Yt with T observations, {Yt = (yt,i,j )n×m : t =1, . . . , T }.
1

We proposed an MAR(p) model with rank-R bilinear forms focusing on the dynamic informationin neighboring columns and rows of Yt ’s to achieve dimension reduction (Hsu, Huang and Tsay, 2021;
Hsu, Huang, Tsay and Kao, 2023):
Yt =

p X
R
X

0
σj,r Aj,r Yt−j Bj,r+ ηt ,(1)

j=1 r=1

which forms a specific case of an vector AR model:yt ≡ vec (Yt ) =

p
X

R
X

j=1

r=1

!
σj,r Bj,r ⊗ Aj,ryt−j + vec(ηt ),(2)

with structured AR coefficient matrices satisfying
Φj =

R
X

σj,r (Bj,r ⊗ Aj,r ) ,j = 1, . . . , p.

(3)

r=1

For model identifiability, we set0
tr(A0j,r Aj,s ) = tr(Bj,rBj,s ) = δrs ,σj,1 ≥ · · · ≥ σj,R > 0;

1 ≤ r, s ≤ R, j = 1, . . . , p,(4)

where δrs = 1 if r = s; and 0 otherwise. The AR coefficient matrix Φj in (3) has a close connection to
the singular value decomposition (SVD) of a re-shaped matrix with rank R.
Regarding the specification of the covariance Ψ for the innovation ηt , we adopt the fixed-rankcovariance model (Cressie and Johannesson, 2008):
Ψ ≡ var(vec(ηt )) = F M F 0 + ση2 I,(5)where F is a pre-specified nm × K matrix of basis functions, M and ση2 ≥ 0 are unknown parameters.
We recommend using the multi-resolution spline basis functions of Tzeng and Huang (2018) for F which
are flexible to characterize features up to a chosen scale.

2

Maximum Likelihood Estimation

Denote Φ ≡ [Φ1 , . . . , Φp ], Yt1 :t2 ≡ [yt1 , · · · yt2 ], and let Yp:(T −1)
yp yp+1 · · · yT −1 ..
..
..  Y(p−1):(T −2) 
..



.
.
. 
Y∗≡ .
=.
..
y2 y3 · · · yT −p+1  .
y1 y2 · · · yT −p
Y1:(T −p)
The Gaussian log-likelihood of (Φ, Ψ), conditional on Y1:p , can be expressed as
0
i
T −p
1 h
`(Φ, Ψ) = −log(det(Ψ)) − tr Y(p+1):T − ΦY ∗ Ψ−1 Y(p+1):T − ΦY ∗ .
2
2
2

(6)

Thus, the constrained ML estimators (MLEs) of (Φ, Ψ) (conditional on Y1:p ) are given by


h
 −1
i
1
∗ 0
∗
(Φ̂, Ψ̂) ≡ arg min log(det(Ψ)) +,tr Y(p+1):T − ΦYΨY(p+1):T − ΦY
T −p
Φ,Ψ

(7)

subject to the constraints (3), (4), (5).

2.1ADMM Algorithm for Obtaining MLEs

The constrained parameter estimates Φ̂ and Ψ̂ do not have closed-form expressions, We develop analternating direction method of multipliers (ADMM) algorithm (Boyd, et al., 2011) to solve the resulting
constrained optimization problem. First, we augment an auxiliary (nm) × (pnm) parameter matrix,C ≡ [C1 , . . . , Cp ] with
Cj =

R
X

σj,r (Bj,r ⊗ Aj,r ) ;j = 1, . . . , p,

(8)

r=1

to alleviate the direct handling of the rank-R structure in Φ. This transfer of the rank-R constraint
from Φ to C simplifies the overall optimization process. By introducing C, the augmented Lagrangianof (7) with the constraint C = Φ follows
n
0
o
1
tr Y(p+1):T − ΦY ∗ Ψ−1 Y(p+1):T − ΦY ∗
T −p
 0
 ρ
+ tr Λ (Φ − C) + kΦ − Ck2F ,2Lρ (Φ, C, Ψ, Λ) = log(det(Ψ)) +(9)

where Λ ≡ [Λ1 , . . . , Λp ] is the matrix of Lagrangian multipliers, ρ > 0 is the augmented Lagrangianparameter. Starting with some initial estimates C (0) , Ψ(0) , Λ(0) of C, Ψ, Λ, the proposed ADMM
algorithm iteratively cycles through four steps: Φ-step, C-step, Ψ-step, and Λ-step, until convergence.
The algorithm is computationally efficient due to the presence of simple closed-form expressions in all
four steps, as detailed in Hsu, Huang, Tsay and Kao (2023). The algorithm’s pseudo-code is outlinedin Algorithm 1.

2.2

Model Selection

The proposed MAR model consists of four structure parameters to select, including the autoregressiveorder p, the reshaped matrix rank R for {Aj,r } and {Bj,r }, the dimension K for the spatial covariancematrix M , and the bandwidth L of {Aj,r } and {Bj,r }. Note that, the rank R of the reshaped matrix
is upper bounded by Rmax ≡ min{n2 , m2 } corresponding to the maximum rank of SVD representationof R(Φj ). We consider Akaike’s information criterion (AIC, Akaike, 1973) and Bayesian information
criterion (BIC, Schwarz, 1978) to select these structure parameters (p, R, K, L):22 
df(Φ) + df(Ψ) ,AIC = −`(Φ̂, Ψ̂) +
T −p
T −p
2
log(T − p) 
df(Φ) + df(Ψ) ,BIC = −`(Φ̂, Ψ̂) +
T −p
T −p
3

Algorithm 1 An ADMM Algorithm for the MLEs
Input:
initial estimates: C (0) , Ψ(0) , Λ(0) ;learning rate ρ;
Output:set k ← 0;
repeatΦ-step: Φ(k+1) = arg min Lρ Φ, C (k) , Ψ(k) , Λ(k) ;ΦC-step: C (k+1) = arg min Lρ Φ(k+1) , C, Ψ(k) , Λ(k) , subject to (4) and (8);C
Ψ-step: Ψ(k+1) = arg min Lρ Φ(k+1) , C (k+1) , Ψ, Λ(k) , subject to (5);Ψ
Λ-step: Λ(k+1) = Λ(k) + ρ Φ(k+1) − C (k+1) ;k ← k + 1;
until a stopping criterion is satisfied.
Table 1: Model selection results for the wind-speed data.
p
1
2
3
4
5
6
7

R
3
3
2
3
2
2
2

K
150
140
140
140
140
140
140

L
3
2
2
1
1
1
2

BIC
-447.13-459.33-462.24-466.02-466.54-469.08-466.86

PMSE
1.556
1.583
1.560
1.546
1.549
1.558
1.534

where df(Φ) and df(Ψ) represent the numbers of free parameters in Φ and Ψ. Under the banded
constraints on (A, B) with a bandwidth L,
df(Φ) = pR[(2n − 1 − L)L + (2m − 1 − L)L + n + m − R], 0 ≤ L ≤ Lmax , 1 ≤ R ≤ Rmax .
Under the Toeplitz constraints on (A, B) with a bandwidth L,
df(Φ) = pR[4L + 1 − (R − 1)],providedR < 4L + 2.

The number of free parameters in Ψ has been shown by Tzeng and Huang (2018) to be:K(K + 1)/2 + 1,
if K < T ,df(Ψ) =
KT + 1 − T (T − 1)/2, if K ≥ T ,where K ≤ nm.

3

Applications

We apply the extended MAR models to the wind-speed data, which consist of 6-hour averaged U-wind(i.e., eastward wind component) values at 289 locations (17x17 regular grid), collected from November
1992 to February 1993 over the region within latitudes 14◦ S-16◦ N and longitudes 145◦ E-175◦ E.
4

Similarly to previous studies (Hsu, Chang and Huang, 2012; Hsu, Huang and Tsay, 2021), we firstremoved the temporal average for each of the 289 locations and applied our method to the mean-adjusted
data {Yt : t = 1, . . . , 480} with Yt being a 17 × 17 matrix.
For fitting and performance evaluation, we split the data into a training sample {Y1 , . . . , Y400 } anda validation sample {Y401 , . . . , Y480 }. We select the best MAR model via BIC among the candidate set:(p, R, K, L) ∈ M ≡ {1, . . . , 7} × {1, . . . , 4} × {10, 20, . . . , 280} × {1, . . . , 16}.
Table 1 reports the performance of our best model compared with the popular models in highdimensional time series. Our MAR model has an overall satisfactory performance among others.

References
Boyd, S., Parikh, N., Chu, E., Peleato, B. and Eckstein, J. (2010). Distributed optimization andstatistical learning via the alternating direction method of multipliers. Foundations and Trends
in Machine Learning, 3(1) 1–122.
Cressie, N. and Johannesson, G. (2008). Fixed rank kriging for very large spatial data sets, Journalof the Royal Statistical Society, Series B, 70, 209–226.
Hsu, N.-J., Huang, H.-C., and Tsay, R. S. (2021). Matrix autoregressive spatio-temporal models,
Journal of Computational and Graphical Statistics, 30, 1143–1155.
Hsu, N-J, Huang, H-C, and Tsay, R. S., Kao, T-C (2023+) Rank-R matrix autoregressive models formodeling spatio-temporal data, forthcoming in Statistics and Its Interface.
Tzeng, S. and Huang, H.-C. (2018). Resolution adaptive fixed rank kriging, Technometrics, 60, 198–
208.
Van Loan, C. and Pitsianis, N. (1993). Approximation with Kronecker products. In Linear Algebra forLarge Scale and Real Time Applications (M. S. Moonen and G. H. Golub eds.), Kluwer Academic
Publishers, Dordrecht, 293–314.

5

Stable P -value Assignment in High-Dimensional Regression via
Data Splitting
Chien-Chung Wanga , Chih-Hao Changb , Hsin-Cheng HuangcabDepartment of Statistics, Colorado State University, USA
Department of Statistics, National Chengchi University, TaiwancInstitute of Statistical Science, Academia Sinica, Taiwan

Abstract
In the rapidly evolving landscape of statistical research, high-dimensional regression models have emerged as a critical tool for tackling complex data structures. These models, whichfeature a large number of predictors relative to the sample size, present unique challenges thatare not adequately addressed by traditional statistical methods. Conventional approaches,such as the ordinary least squares (OLS) estimator, often fall short of providing reliable
estimates and inferences in high-dimensional settings. This limitation is particularly evident
in the realms of variable selection, confidence interval construction, and hypothesis testingtargeting individual regression coefficients.
A plethora of specialized model selection algorithms have been developed to navigate thecomplexities of high-dimensional data. These include Lasso, SCAD, Adaptive Lasso, SIS,
MCP, and OGA, among others. Each of these algorithms brings its own set of advantages
and limitations to the table, and they have been widely adopted in various applicationsof high-dimensional regression analysis. However, a common challenge that persists acrossthese methods is the difficulty in conducting rigorous statistical inference on the selectedmodels. Traditional inferential procedures often fail to account for the intricacies involved
in high-dimensional model selection, leading to biased or inconsistent results.
In this work, we introduce a series of pivotal modifications to the existing multiplesplitting paradigm, aiming to improve the robustness and reliability of variable selection in
high-dimensional regression models. Our modifications are designed to mitigate the conser-

1vatism that is inherent in the original multiple-splitting approach. We redefine the p-value
assignment mechanism for individual data splits, incorporating a more nuanced statisticalframework that better captures the underlying data distribution. Furthermore, we integrate
the stability-selection approach into our methodology, which allows for a more rigorous control of the false discovery rate.
To further enhance the overall performance of the variable selection process, we introduce the concept of multiple dependent p-value aggregation. This innovative framework
capitalizes on the underlying correlation structure among the predictors, thereby offeringa more comprehensive and accurate variable selection strategy. Unlike existing methods,
our approach does not focus on computational efficiency but prioritizes the robustness andreliability of the variable selection process.
To validate our methodology, we conduct extensive simulation studies encompassing various scenarios, including varying sample sizes, predictor correlations, and signal-to-noise ratios. These simulations confirm the robustness of our approach in different high-dimensionalsettings. While we have not yet applied our methodology to real-world data sets in fields
such as finance and social sciences, the empirical results suggest that it holds promise forsuch applications.
In summary, our work offers an advancement in p-value calculations for high-dimensional
regression models. By redefining the p-value assignment mechanism and introducing multiple dependent p-value aggregation, we provide a more rigorous and reliable framework forstatistical inference. This focus on p-values ensures that our methodology not only excels invariable selection but also provides more accurate and robust hypothesis tests. We believethese advancements will be a cornerstone for future research in high-dimensional statisticalinference.
Keywords: Exchangeable p-values; Family-wise error rate; False discovery rate; Multiple
testing; Post-selection inference; Selection consistency; Variable selection2Modeling Joint Cylindrical Distributions
and Related Markov Processes
Toshihiro Abe1 , Tomoaki Imoto2 ,Takayuki Shiohama3 , Yoichi Miyata41

Hosei University, 2 University of Shizuoka,
3
Nanzan University 4 Takasaki City University of Economics

Abstract
Joint cylindrical distributions are key probability distributions that make possible a multivariate regression analysis as well as Markov models on a cylinder. In this study, somejoint cylindrical distributions are proposed, and their statistical properties together withalgorithms for random number generation are investigated. For proposed joint cylindrical
distributions, the marginal distributions for various combinations of linear and circular variables are obtained. These marginal distributions are applied to the Markov process on acylinder. The maximum likelihood estimation for unknown model parameters is investigated. To illustrate the applicability of the proposed models, time series analysis using windspeed and direction data is investigated.
keywords : circular statistics; cylindrical data; Markov process; maximum likelihood estimation.

1

Introduction

Bivariate data which has a form of the combination of angular and linear variables are calledcylindrical data. Examples of cylindrical data can be found in various ﬁelds in the environment
and natural science, engineering, information science, and so forth. Because a cylindrical data
analysis may provide clear insights into the probabilistic structures and relationships between
the angular and linear variables, many common and conventional statistical tools are applied tothe data analysis.
Fitting cylindrical data using cylindrical distribution has a long history. Example of earlyworks includes the distribution proposed by [5], where they considered the exponential distribution on a linear variable and von Mises distribution on a circular variable. The normal
distribution on a linear part and von Mises distribution on a circular part are proposed by [7],where the mean and variance of the normal distributions are the functions of circular variables.
The drawbacks of this model were found to be diﬃcult to express the marginal distribution ofa linear variable. In [2], a Weibull von Mises distribution on a cylinder was considered andillustrated the model’s tractability for real data analysis. This distribution was extended todescribe heavy tails on linear variables by [4].
Time series analysis for cylindrical data is of special importance in the prediction of windspeed and directions in wind ﬁrm management with accurately predicting wind energy. Hidden
Markov Models based on the Weibull von Mises distribution proposed by [2] are considered by[6], and they found that their models can explain parsimonious accommodations of circular-linear
correlation, multimodality, skewness, and temporal autocorrelation of the cylindrical data.

Joint cylindrical distributions are beneﬁcial as they provide a regression analysis for both
response and explanatory variables that take values on a cylinder. For time series modeling, thejoint cylindrical distribution is needed to formulate a Markov model on a cylinder. Despite thenecessity for such needs in statistical modeling, there is little known about the joint distributionof a cylinder. Recently, joint cylindrical distribution is proposed by [1]. However, the proposed
distribution has several shortcomings to ﬁt the models in real data applications. The idea
behind this study is a natural extension of the modeling that is used in the study [3] and [1].
The possible way of constructing a joint cylindrical distribution is to apply an approach adopted
by [8] and [5]. However, this modeling needs to handle distribution functions for each marginaldensity which makes it diﬃcult to formulate various combinations of joint marginal distributions.
Similar diﬃculties arise when we apply copula-based multivariate modeling.
The joint cylindrical distributions proposed in this study and their applications to regression
analysis and time series analysis are classiﬁed into a distributional regression approach. Thesestatistical methods may provide more ﬂexible modeling than nonparametric regression models
as they provide whole probabilistic information on the conditional distribution of the response
variable given explanatory variable, including quantiles and median of the predicted response
variable.
Because of its mathematical tractability and a good ﬁt to real-world data, the Markov
property of a conditional cylindrical distribution is a basic concept in state-space modeling andis widely applied in various ﬁelds. In this study, we elucidate how to construct joint cylindrical
distributions that are tractable for data analysis and deduce conditional distribution which isapplicable for Markov transition density for time series analysis. We show that the resulting
Markov processes are found to have a non-linear ﬁrst-order autoregressive structure. In addition,
conditional circular autocorrelation functions are obtained for the proposed processes.

References[1] Abe, T., Imoto, T., Shiohama, T. & Miyata, Y. : On some ﬂexible models for circular,
toroidal, and cylindrical data. Directional Statistics for Innovative Applications, Springer,
pp. 229–243. (2022)[2] Abe, T. & Ley, C. : A tractable, parsimonious and ﬂexible model for cylindrical data, withapplications. Econometrics and Statistics, 4, pp. 91-104. (2017)[3] Imoto, T. & Abe, T. : The simple construction of a toroidal distribution from independentcircular distributions. Journal of Multivariate Analysis, 186, 104799 (2021)[4] Imoto, T., Shimizu, K., & Abe, T. : A cylindrical distribution with heavy-tailed linear
part. Japanese Journal of Statistics and Data Science, 2 (1), pp. 129–154. (2019)[5] Johnson, R. A. & Wehrly, T. E. : Some angular-linear distributions and related regression
models. Journal of the American Statistical Association, 73 (363), pp. 602–606. (1978)[6] Lagona, F., Picone, M. & Maruotti, A. : A hidden Markov model for the analysis ofcylindrical time series. Environmetrics, 26 (8), pp. 534–544. (2015)[7] Mardia, K. V. and Sutton, T. W. : A model for cylindrical variables with applications.
Journal of the Royal Statistical Society: Series B, 40 (2), pp. 229–233. (1978)[8] Wehrly, T. E. & Johnson, R. A. : Bivariate models for dependence of angular observations
and a related Markov process. Biometrika, 67 (1), pp. 255–256. (1980)
2

Test for the existence of the residual spectrum withapplication to brain functional connectivity detection
Yuichi Goto
Xuze Zhang
Benjamin Kedem
Shuo Chen(Kyushu University)∗1(University of Maryland)(University of Maryland)(University of Maryland)

本講演では, 非線形構造を考慮した残差スペクトルという時系列解析の周波数領域における指標を導入し, 残差スペクトルが存在するかどうかの検定問題に対する検定手法を提案した. コヒーレンスは, ２つの時系列間の類似度の指標であり, ピアソンの相関係数の時系列への拡張となっている. しかし, コヒーレンスでは時系列間の線形関係しか捉えることができないため, 非線形関係を捉えられる指標が必要とされている. 応用
例として, 脳障害と深く関連している脳の機能的結合の検出が挙げられ, コヒーレンス
では検出できなかった機能的結合が残差スペクトルでは検出が可能であることを示した. この講演は, Goto et al. (2023) に基づいたものである.

1. モデルと主定理
本講演では, 次のモデルを考えた;

X0 (t) := ζ +

K
∞
X
X

bi (ki )Xi (t − ki ) + ϵ(t).

(1)

i=1 ki =−∞ここで, X0 (t) は応答変数, ζ は切片, X(t) := (X1 (t), . . . , XK (t))⊤ は共変量であり, 自己共分散行列 ΓX (u) := EX(t)X(t − u) を持つ. ただし, Xi (t) := Xi′ − EXi′ であり,′(t))⊤ は, 任意の s ∈ N に対して, s 次定常過程. また, ϵ(t) は任意の(X1′ (t), . . . , XK
モーメントを持つ平均ゼロである i.i.d. 誤差過程で {X; t ∈ Z} と独立, 係数 bi (k) はP∞
2
k=−∞ k |bi (k)| < ∞ を満たす.
次の条件を仮定する:
仮定 1. 任意の２以上の整数 ℓ と (i1 , . . . , iℓ ) ∈ {1, . . . , K}ℓ に対して,
!
∞
ℓ
X
X
2
1+
|sj | |Cum{Xi1 (0), Xi2 (s2 ), . . . , Xiℓ (sℓ )}| < ∞s2 ,...,sℓ =−∞

j=2

を満たす.
講演では, 適当なモーメント条件を満たす geometrically α-mixing s 次定常過程は, 仮
定 1 を満たすということを示した.
仮定 1 の下で, (X0 (t), X1 (t), . . . , XK (t)) に対するスペクトル行列 f (λ) := (fij (λ))i,j=0,...,Kが存在し, 二回連続微分可能である. モデル (1) の直交分解の結果, 次で定義される残差This research was supported by JSPS Grant-in-Aid for Early-Career Scientists JP23K16851 (Y.G.)
and NIH grant DP1 DA048968 (S.C.). We thank Brin Mathematics Research Center for supporting
the ﬁrst author’s (Y.G.) research stay at the University of Maryland. Part of this research was
conducted while the ﬁrst author (Y.G.) was visiting the Brin Mathematics Research Center at the
University of Maryland.
∗1e-mail: yuichi.goto@math.kyushu-u.ac.jp

スペクトルが自然に導出されることを明らかにした:2
Ajj eiλ det (fj (λ))|f10 (λ)|2
fG1 G1 (λ) =and fGj Gj (λ) =f11 (λ)det (fj−1 (λ))

for j = 2, . . . , K.

†
ここで, fj (λ) := (fij (λ))i,j=1,...,j と fa,b:= (f1b (λ), . . . , fab (λ))⊤ に対して,


Pj−1
‡
−
det
f
(λ)
fi0 (λ) + det (fj−1 (λ)) fj0 (λ)

i,j
i=1
Ajj eiλ :=,det (fj (λ))


‡
†
†
†
†
†
fi,j
(λ) := fj−1,1(λ), . . . , fj−1,i−1(λ), fj−1,j(λ), fj−1,i+1(λ), . . . , fj−1,j−1(λ) .

上で定義した AKK eiλ の分子を ΦK (f (λ)) と書くとする. 次の仮説検定問題を考
える:

Z π

|ΦK (f (λ))|2 = 0, v.s. K0 :H0 が成立しない.

H0 :
−π

これは, モデル (1) に {XK (t)} が含まれているかどうかに対応する. 検定統計量を

n
Tn := √
Mn

Z π


2
ΦK fˆ(λ) dλ − µ̂n,K−πと定義する. ここで, fˆ(λ) は f (λ) のカーネル推定量, µ̂n,K はバイアス調整項で

Z



π

√
µ̂n,K := M n ηω,2



fˆ(λ)





 ∂ΦK
tr 
∂Z ⊤
−π

fˆ(λ)
Z=fˆ(λ)

∂ΦK

fˆ(λ)






fˆ(λ) dλ,∂Z
Z=fˆ(λ)R∞ηω,2 := −∞ ω 2 (x)dx < ∞, ω(x) はカーネル推定量に含まれるラグウィンドウ. 仮定 12とその他の正則条件, 帰無仮説 H0 の下で, Tn は平均ゼロで分散 σKの正規分布に分布収
R∞ 4束する. ここで, ηω,4 := −∞ ω (x)dx < ∞,Z

2σK:= 4πηω,4πtr∂ΦK (f (λ))∂ΦK (f (λ))
f (λ)
f (λ)
⊤
∂Z
∂Z
Z=f (λ)Z=f (λ)

!2

−π

+ tr

∂ΦK (f (λ))∂ΦK (f (λ))
f (λ)
f (λ)
⊤
∂Z
∂Z ⊤
Z=f (λ)Z=f (λ)!2
dλ.

したがって, 有意水準 α に対して, Tn /σ̂K ≥ zα を満たすとき帰無仮説 H0 を棄却する
22
検定は漸近サイズ α である. ここで, zα は標準正規分布の上側 α%点, σ̂K
は σK
の一致
推定量.
さらに講演では提案検定手法が 一致性を持つことも証明した. 講演では, 数値実験の結果と脳の機能的結合への応用例について紹介した.

References
Goto, Y., Zhang, X., Kedem, B., and Chen, S. (2023). Residual Spectrum Applied inBrain Functional Connectivity. ArXiv E-prints, Available at arXiv:2305.19461.

Prediction-based statistical inference formultiple time series
Yan Liu 1
Abstract
A new minimum contrast estimator is proposed in this work for multivariate time series in thefrequency domain. It is based on the prediction problem of time series. The properties of the
proposed contrast estimation function are explained in details. We also derive the asymptotic
normality of the proposed estimator and compare the asymptotic variance with the existing results. The asymptotic eﬃciency of the proposed minimum contrast estimation is also considered.
The theoretical results are illustrated by some numerical simulations.
keywords: Minimum contrast estimation, Prediction problem, Spectral analysis, Multiple timeseries, Spectral density matrix, Asymptotic eﬃciency
1. Prediction-based contrast function
Let {z(t)} be a second-order stationary processz(t) =∞∑

G(j)e(t − j),t ∈ Z.

(1.1)

j=0

Assuming∑∞

j=0 tr G(j)KG(j)⊤f (ω) =

< ∞, the process has a spectral density matrix
1k(ω)Kk(ω)∗ ,2π

ω ∈ Λ := (−π, π],∑∞
where k(ω) = j=0 G(j)eijω .
We propose a new contrast function D for multivariate time series analysis as follows:∫ πD(fθ , g) =a(θ) tr[|fθ (ω)|p g(ω)] dω,
(1.2)
−π

√
where |fθ (ω)| := fθ (ω)∗ fθ (ω) is the square root of the matrix fθ (ω)∗ fθ (ω), and()−p/(p+1)

 ∫π
p+1
tr[|f
(ω)|
]
dω
,
p ̸= −1,θ−πa(θ) =1,p = −1.
Lemma 1.1. Let p < 0. Suppose the matrices fθ and g are both Hermitian and semi-deﬁnite
in any ω ∈ Λ. The following inequality holds:)1/(p+1)(∫ πD(fθ , g) ≥tr[|g(ω)|p+1 ] dω.

−π

1 This work was supported by JSPS Grant-in-Aid for Scientiﬁc Research (C) 23K11018.

1

(1.3)

Yan Liu/ Prediction-based statistical inference2

2. Statistical inference
In this section, we discuss the parameter estimation problem based on the new contrast function(1.2). The following assumptions are imposed.
Assumption 2.1.
(i) The parameter space Θ is a compact subset of Rd .
(ii) If θ1 ̸= θ2 , then fθ1 ̸= fθ2 on a set of positive Lebesgue measure.
(iii) The spectral density matrix fθ (ω) is three times continuously diﬀerentiable with respect∂2
to θ and the second derivative ∂θ∂θ⊤ fθ (ω) is continuous in ω.
Lemma 2.2.
Under Assumption 2.1, when p < 0, θ0 minimizes the disparity D(fθ , fθ0 ) if θ0 ∈ Θ.
Theorem 2.3. Under regularity conditions, if θ0 ∈ Θ ⊂ R, then the following holds:
(a)
(b)

p

θ̂n −→ θ0 .
√dn(θ̂n − θ0 ) −→ N (0, H(θ0 )−1 V (θ0 )H(θ0 )−1 ).

Assumption 2.4.
• The Gaussian-Fisher information matrix F (θ0 ) is invertible.
• The eigenvectors of fθ (λ) do not depend on the parameter θ.
Theoretically, it can be shown under assumptions that
F (θ0 )−1 ≤ H(θ0 )−1 Ṽ (θ0 )H(θ0 )−1 ,where the equality holds then α = −1, or the spectral density matrices does not depend on ω.
3. Numerical simulations
In this section, we present the numerical performance of our minimum contrast estimator for
multivariate time series. The data are generated from the following models:Model (i)

(
Xt =

)
0.2 θ1
Xt−1 + ϵt ,0 0.4

ϵt ∼ N (0, I2×2 );Model (ii)
(
Xt =

)
(
)
0.2 θ200
Xt−1 +
Xt−1 + ϵt ,0 0.4
θ2 /3 0ϵt ∼ N (0, I2×2 ),where I2×2 denotes the identity matrix of size 2.
The ﬁnite sample performance of the new minimum contrast estimators is shown by the Q-Q
plots, and a thorough investigation based on diﬀerent parameters θ1 = 0, 0.2, 0.4, 0.6, 0.8, 1 andθ2 = 0, 0.2, 0.4, 0.6, 0.8, 1. The size of the data is n = 100. We conclude that the estimators whenp = −2, −3, −4 concentrate around 0 more than the Whittle estimator (when p = −1).
References
Liu, Y. (2023). A Minimum Contrast Estimation for Spectral Densities of Multivariate Time
Series. In Research Papers in Statistical Inference for Time Series and Related Models 325–
342. Springer.

Weak convergence of the partial sum of I(d) process to a
fractional Brownian motion in ﬁnite interval representation
Junichi Hirukawa and Kou Fujimori
Niigata University and Shinshu University
ABSTRACTAn integral transformation which changes a fractional Brownian motion to a process with independent increments
has been given. A representation of a fractional Brownian motion through a standard Brownian motion on a ﬁnite
interval has also been given. On the other hand, it is known that the partial sum of the discrete time fractionallyintegrated process (I(d) process) weakly converges to a fractional Brownian motion in inﬁnite interval representation.
In this talk we derive the weak convergence of the partial sum of I(d) process to a fractional Brownian motion in ﬁnite
interval representation.

1 Introduction
Stochastic analysis for FBM has been developed by Decreusefond and Üstünel (1997) using Malliavin calculus. Norros
et al. (1999) showed that many basic results can be obtained more directly with rather elementary arguments andcomputations. Norros et al. (1999) considered a normalized fractional Brownian motion (FBM) (Zt )g≥0 with selfsimilarity parameter H ∈ (0, 1). Mandelbrot and Van Ness (1968) deﬁend the process more constructively as theintegral
!
Z t
Z t n
o
(t − u)H−1/2 dWu +(t − u)H−1/2 − (s − u)H−1/2 dWu ,Zt − Z s = c H
−∞

s

 
where Wt is the standard Brownian motion. The normalization E Z12 = 1 is achieved with the choice



1/2


2HΓ 32 − H

 ,

cH =  Γ H + 12 Γ (2 − 2H)
where Γ (·) denotes the Gamma function.

1.1 The fundamental martingale M
Norros et al. (1999) considered the following process. Let w (t, s) be the function


1/2−H


(t − s)1/2−H , for s ∈ (0, t) ,c1 sw (t, s) = 0,for s < (0, t) ,
where
!)−1
(
1
1
− H, H +c1 = 2HB
2
2

1

and B is the beta funtion
B (u, v) =

Γ (u) Γ (v).
Γ (u + v)

Then, the centered Gaussian process
Mt =

Z t
w (t, s) dZ s0

has independent increments and variance function 
E Mt2 = c22 t2−2H ,
where
c2 =

cH
.
2H (2 − 2H)1/2

In particular, M is a martingale.

2 Weak convergence of I (d) process
Now, we obtain the following functional central limit result for I (d) process, if d > 0,
[nt]−1[nt][nt]−1

X 
X


1 e
1
1
X
1/2
1/2 
v
DW
+
θ
v
DW s
Z
=


[nt]
s
u,u+1−s


s−1
s−1

σnd+1/2
σnd+1/2 s=1
σnd+1/2 s=1  u=s
)
(Z t
Z t
Z t
1
(u − s)d−1 ud du dW (s) :=⇒s−ddZ (s) = Z (t) .
Γ (d) 0s
0References
Mandelbrot B. B. and van Ness J. W. (1968). Fractional Brownian motions, fractional noises and applications SIAM
Review, 10, 422–437.
Norros, I. and Valkeila, E. and Virtamo, J. (1999). An elementary approach to a Girsanov formula and other analyticalresults on fractional Brownian motions Bernoulli 4, 571–587.
Pipiras, V. and Taqqu, Murad S. (2001). Are classes of deterministic integrands for fractional Brownian motion on
an interval complete? Bernoulli 7, 873–897.
Tanaka, K. (2013). Distributions of the maximum likelihood and minimum contrast estimators associated with thefractional Ornstein-Uhlenbeck process Stat. Inference Stoch. Process. 16, 173–192.

2