International Symposium on Theories and Methodologies
for Large Complex Data
November 21-23, 2019
Venue:
Conference Room 406, Tsukuba International Congress Center
2-20-3 Takezono, Tsukuba, Ibaraki 305-0032, Japan
Organizers:
Makoto Aoshima (University of Tsukuba)
Mika Sato-Ilic(University of Tsukuba)

Kazuyoshi Yata

(University of Tsukuba)

Aki Ishii(Tokyo University of Science)

Supported by
Grant-in-Aid for Scientific Research (A) 15H01678 (Project Period: 2015-2019)“Theories and methodologies for large complex data”(Principal Investigator: Makoto Aoshima)
Grant-in-Aid for Challenging Research (Exploratory) 19K22837 (Project Period: 2019-2021)“Tackling individualized modeling with ultra-high dimensional data”(Principal Investigator: Makoto Aoshima)

Program
November 21 (Thursday)
14:00∼14:10

Opening

14:10∼14:50 Aki Ishii∗,a , Kazuyoshi Yatab and Makoto Aoshimabab(Department of Information Sciences, Tokyo University of Science)(Institute of Mathematics, University of Tsukuba)

Tests for high-dimensiomal covariance structures under the SSE model
15:00∼15:40 Takahiro Nishiyama∗,a , Masashi Hyodob and Tatjana Pavlenkoca(Department of Business Administration, Senshu University)b(Department of Mathematical Sciences, Osaka Prefecture University)c(Department of Mathematics, KTH Royal Institute of Technology)

On error bounds for high-dimensional asymptotic distribution of L2 -type test statistic
15:55∼16:35 Hiroumi Misaki(Faculty of Engineering, Information and Systems, University of Tsukuba)

Financial risk management with high-frequency data(∗ Speaker)

16:45∼17:25 Junichi Hirukawa∗,a and Kou Fujimoribab(Faculty of Science, Niigata University)(School of Fundamental Science and Engineering, Waseda University)

Weak convergence of the partial sum of I(d) process to a fractional Brownian motion
in ﬁnite interval representation

November 22 (Friday)
9:20∼10:00

Kengo Kamatani (Graduate School of Engineering Science, Osaka University, and JST CREST)

High-dimensional analysis of the piecewise deterministic Markov process
for Bayesian inference
10:10∼10:50 Shogo Kato∗,a and Peter McCullaghbab(The Institute of Statistical Mathematics)(Department of Statistics, University of Chicago)

A Cauchy family derived by the Möbius transformations of the sphere
11:00∼17:35

Special Invited and Keynote Sessions

18:30∼

DinnerNovember 23 (Saturday)
9:20∼10:00 Shota Katayama(Faculty of Economics, Keio University)

Direct estimation of conditional averaging treatment eﬀect in high dimensions
10:10∼10:50 Kei Hirosea,∗ and Hiroki Masudabab(Institute of Mathematics for Industry, Kyushu University)(Faculty of Mathematics, Kyushu University)

Statistical modeling for electricity load forecasting
11:00∼11:40 Takuma Bandoa , Tomonari Sei∗,a and Kazuyoshi Yatabab(Graduate School of Information Science and Technology, University of Tokyo)(Institute of Mathematics, University of Tsukuba)

Consistency of the objective general index in high dimensional settings
11:40∼ 11:50

Closing(∗ Speaker)

Special Invited Session
11:00∼11:50

Data beyond the euclidean space

Speaker: Jörn Schulz(Department of Electrical engineering and Computer science, University of Stavanger)
Chair: Shogo Kato
11:50∼13:15(The Institute of Statistical Mathematics)

Lunch

13:15∼14:05 Change points detection and identiﬁcation for high dimensional dependent data

Speaker: Ping-Shou Zhong(Department of Mathematics, Statistics, and Computer Science, University of Illinois at Chicago)
Chair: Fumiya Akashi (Graduate School of Economics, University of Tokyo)

14:15∼15:05

Towards a sparse, scalable, and stably positive deﬁnite (inverse) covarianceestimator

Speaker: Joong-Ho (Johann) Won(Department of Statistics, Seoul National University)
Chair: Shota Katayama (Faculty of Economics, Keio University)

Keynote Session
15:20∼16:20

A two-stage dimension reduction method and its applications
on highly contaminated image sets

Speaker: I-Ping Tu(Institute of Statistical Science, Academia Sinica)
Discussion Leader: Yuan-Tsung Chang

16:35∼17:35(Department of Social Information, Mejiro University)

Sample covariance matrices from “bad populations”

Speaker: Jeﬀ Yao(Department of Statistics and Actuarial Science, The University of Hong Kong )Discussion Leader: Kazuyoshi Yata (Institute of Mathematics, University of Tsukuba)Tests for high-dimensiomal covariance structures
under the SSE model
Aki Ishiia , Kazuyoshi Yatab and Makoto Aoshimaba

Department of Information Sciences, Tokyo University of SciencebInstitute of Mathematics, University of Tsukuba

1 Introduction
In this talk, we consider testing the high-dimensional intraclass covariance matrix. We produce a new test
statistic for each covariance structure by using the extended cross-data-matrix (ECDM) methodology. We
show that the test statistic is an unbiased estimator of its test parameter. We prove that the test statistic has
a consistency property and establishes the asymptotic normality. We propose a new test procedure for thehigh-dimensional intraclass covariance matrix and evaluate its asymptotic size and power theoretically.
Suppose we take samples, xj = (x1j , . . . , xpj )T , j = 1, . . . , n, of size n (≥ 4), which are independent and identically distributed (i.i.d.) as a p (≥ 2)-variate distribution. We assume that xj has an
unknown mean vector µ and unknown (positive-semidefinite) covariance matrix Σ. Let σ = tr(Σ)/p.
Let σij be the (i, j) element of Σ for i, j = 1, . . . , p. We assume that σjj ∈ (0, ∞) as p → ∞for all j. For a function, f (·), “f (p) ∈ (0, ∞) as p → ∞” implies that lim∑
inf p→∞ f (p) > 0 andlim supp→∞ f (p) < ∞. Then, it holds that σ ∈ (0, ∞) as p → ∞. Let ρ = pi̸=j σij /{σp(p − 1)}.
Note that
1Tp Σ1p= σ{1 + ρ(p − 1)}
p

(1.1)

and ρ ∈ [−(p − 1)−1 , 1], where 1p = (1, . . . , 1)T . We denote the identity matrix of dimension p by I p .
In this paper, we consider testing
H0 : Σ = Σ∗vs. H1 : Σ ̸= Σ∗ ,(1.2)

where Σ∗ is a candidate (positive-semidefinite) covariance matrix. For Σ∗ we consider the following covariance structures: (i) identity matrix, (ii) scaled identity matrix, (iii) diagonal matrix, and (iv) intraclasscovariance matrix. Let
ΣS = σI p , ΣD = diag(σ11 , . . . , σpp ) and ΣIC = σ{(1 − ρ)I p + ρ1p 1Tp }.
Let
Σ0 = Σ − Σ∗ and ∆ = ∥Σ0 ∥2F = tr(Σ20 ),where ∥ · ∥F is the Frobenius norm. Note that ∆ = 0 under H0 and ∆ > 0 under H1 . We regard ∆ as atest parameter and construct a test procedure for (1.2) by using an estimator of ∆.
In the current paper, we take a different nonparametric approach and produce a new test statisticfor (1.2). We utilize the extended cross-data-matrix (ECDM) method developed by Yata and Aoshima[2] which is an extension of the cross-data-matrix methodology created by Yata and Aoshima [1]. The
ECDM method is a nonparametric method to produce an unbiased estimator for a function of Σ at a low
computational cost even for ultra high-dimensional data.

2 Unbiased estimator of ∆∑Let Aj be a p × p known idempotent matrix with rank rj (≥ 1) for j = 1, . . . , q, such that qj=1 rj = p∑
and qj=1 Aj = I p , where r1 ≤ · · · ≤ rq when q ≥ 2. Note that tr(Aj ) = rj , A2j = Aj andAj Aj ′ = O for all j (̸= j ′ ). Let κj (≥ 0) be an unknown scalar such that tr(ΣAj ) = rj κj for all j.
Hereafter, we assume that Σ∗ has the following structure:

Note that tr(Σ2∗ ) =as follows:∑qΣ∗ = κ1 A1 + · · · + κq Aq .

(2.1)

2
2
2
2
2
j=1 rj κj and ∆ = tr(Σ ) − tr(Σ∗ ), so that tr(Σ ) ≥ tr(Σ∗ ). One can summarize(I) A1 = I p , κ1 = σ, r1 = p and q = 1 when Σ∗ = ΣS ;(II) Aj = diag(0, . . . , 0, 1, 0, . . . , 0) whose j-th diagonal element is 1, κj = σjj , rj = 1 for all j andq = p when Σ∗ = ΣD ;(III) A1 = 1p 1Tp /p, A2 = I p − 1p 1Tp /p, κ1 = σ{1 + (p − 1)ρ}, κ2 = σ(1 − ρ), r1 = 1, r2 = p − 1
and q = 2 when Σ∗ = ΣIC .
We note that tr(ΣAj )2 /rj = rj κ2j for all j.
We first give an unbiased estimator of ∆ by using the ECDM method. Let n(1) = ⌈n/2⌉ andn(2) = n − n(1) , where ⌈x⌉ denotes the smallest integer ≥ x. Let{{⌊k/2⌋ − n(1) + 1, . . . , ⌊k/2⌋}
if ⌊k/2⌋ ≥ n(1) ,V n(1)(k) ={1, . . . , ⌊k/2⌋} ∪ {⌊k/2⌋ + n(2) + 1, . . . , n} otherwise;{{⌊k/2⌋ + 1, . . . , ⌊k/2⌋ + n(2) }if ⌊k/2⌋ ≤ n(1) ,V n(2)(k) ={1, . . . , ⌊k/2⌋ − n(1) } ∪ {⌊k/2⌋ + 1, . . . , n} otherwise
for k = 3, . . . , 2n−1, where ⌊x⌋ denotes the largest integer ≤ x. Let #S denote the number of elements
in a set S. Note that #V n(l)(k) = n(l) , l = 1, 2, V n(1)(k) ∩ V n(2)(k) = ∅ and V n(1)(k) ∪ V n(2)(k) ={1, . . . , n} for k = 3, . . . , 2n − 1. Also, note that i ∈ V n(1)(i+j) and j ∈ V n(2)(i+j) for i < j (≤ n).
Let
∑
∑
−1
x(1)(k) = n−1
x
=
n
x
and
xj
j
(2)(k)(1)(2)
j∈V n(1)(k)j∈V n(2)(k)for k = 3, . . . , 2n − 1. Let un(l) = n(l) /(n(l) − 1) for l = 1, 2,y ij(1) = xi − x(1)(i+j)

and y ij(2) = xj − x(2)(i+j)

for all i < j. We note that un(l) E(y ij(l) y Tij(l) ) = Σ for l = 1, 2, and y ij(1) and y ij(2) are independent
for all i < j. For example, Yata and Aoshima [2] gave an estimator of tr(Σ2 ) asWn =n
2un(1) un(2) ∑(y Tij(1) y ij(2) )2n(n − 1)

(2.2)

i<j

by the ECDM method. Then, it holds that E(Wn ) = tr(Σ2 ). We also give an unbiased estimator oftr(Σ2∗ ) by using the ECDM method.

References[1] Yata, K. and Aoshima, M. (2010). Effective PCA for high-dimension, low-sample-size data withsingular value decomposition of cross data matrix. Journal of Multivariate Analysis, 101, 2060–
2077.
[2] Yata, K. and Aoshima, M. (2013). Correlation tests for high-dimensional data using extended crossdata-matrix methodology. Journal of Multivariate Analysis, 117, 313–331.

On error bounds for high-dimensional asymptotic distributionof L2 -type test statistic
Takahiro Nishiyamaa , Masashi Hyodob and Tatjana Pavlenkocab

Department of Business Administration, Senshu University
Department of Mathematical Sciences, Osaka Prefecture Universityc
Department of Mathematics, KTH Royal Institute of Technology

This talk was concerned with the canonical testing problem in modern statistical inference, namely the two-sample test for equality of means of independent multivariate populations with very large dimensions. Precisely, let xgk = (xg1k , . . . , xgpk )> , k ∈ {1, . . . , ng },p×p
be ng iid random vectors with xgk » Np (µg , Σg ), g ∈ {1, 2}, where µg ∈ Rp , Σg ∈ R>0,
represent the usual parameters.
We are interested in testing the hypothesis H : kµ1 ¡ µ2 k = 0, vs. A : kµ1 ¡ µ2 k > 0
where k ¢ k denotes L2 -norm. The best-known test procedure which accommodates highdimensional data and allows for Σ1 6= Σ2 is the L2 -type test statistics introduced by Chen
and Qin (2010) (hereafter called for Ch-Q test):
Tn = kx1 ¡ x2 k2 ¡2X

tr(Sg )/ng ,g=1

where n = n1 + n2 , xg and Sg are the sample mean and sample covariance matrix ofgth population. To the unbiasedness property of Tn , Chen and Qin (2010) quantiﬁed its
variance
2X
2tr(Σ2g )4tr(Σ1 Σ2 )2σn = var(Tn ) =+,n (ng ¡ 1)
n1 n2
g=1 g
under H, and show that the distribution of Tn admits a normal limit after appropriate
rescaling. We ﬁrst presented a theoretical analysis of this problem which requires understanding of the rate of convergence of Tn to its normal limit, and then introduced two
new approximations which are more accurate in the asymptotic regime where both n andp tend to inﬁnity.
To provide the intuition behind the approximations which we propose, let ψTen (t)denote the characteristic functions of Ten = Tn /σn . Also, let λr (Λ) be the r-th largest
eigenvalue of the matrix Λ = (n/n1 )Σ1 + (n/n2 )Σ2 and let ∆ = λ1 (Λ)2 /tr(Λ2 ).
Then we obtained the characteristic function of Ten as
ψTen (t) ≈ e¡t /2 + a(it)3 e¡t /2
22

where a = 4b/(3σn3 ) andb=

2X(ng ¡ 2)tr(Σ3g )
g=1

n2g (ng ¡ 1)2+

3tr(Σ21 Σ2 ) 3tr(Σ1 Σ22 )+.
n21 n2n1 n22(1)

By this result, the normal approximation of Ten is immediately achieved, that is FTen (t) =Φ(t) + o(1), as ∆ = 0. On the other hand, inverting the right side of (1) term by term
provided the approximating distribution of Ten of the form F2 (x) = Φ(x) + a(1 ¡ x2 )φ(x),where φ(x) denotes the density of Φ(x).
Another avenue of research delivers a χ2 -approximation of Ten ; this look on the problem
is motivated by the work of Buckley and Eagleson (1988) and Zhang (2005). Let Gd (¢)
denote the cumulative distribution function of χ2 -distribution
with d degrees of freedom.
√
Then we could approximate FTen (x) by F3 (x) = Gd ( 2dx + d).
To establish the rate of convergence, we obtained the explicit bounds for the Kolmogorov distance between the distribution of Ten and its approximations Φ(x), F2 (x), andF3 (x).
Theorem 1. The distribution of Ten satisﬁes the following properties under H:(i) For any n1 , n2 , p, Σ1 , Σ2 such that ∆ < 1/8,
½
¾
3∆
8!1/4
8∆(2 + ω)
sup |FTen (x) ¡ F2 (x)| •
2+
+
.
2πω
8(1 ¡ 8∆)2
9πω 2
x∈R
(ii) For any n1 , n2 , p, Σ1 , Σ2 such that ∆ < 1/8,
½
¾
3∆
8!1/4
{10 + 3(1 ¡ 8∆)¡2 }∆sup |FTen (x) ¡ F3 (x)| •
2+
+
.
2πω
8(1 ¡ 8∆)22πx∈R(iii) For any n1 , n2 , p, Σ1 , Σ2 such that ∆ < 1/6,
(
)
√
√
2∆1/2
6!1/4 2
sup |FTen (x) ¡ Φ(x)| •
3 π+
,
2πω
3(1 ¡ 6∆)3/2
x∈R
where ω is omega constant which is a mathematical constant deﬁned by ω exp(ω) = 1, andω ≈ 0.56714.
Based on Theorem 1, we established convergence rate of the proposed approximations.
Also, we proposed the Ch-Q test with the proposed approximations rests on the adjusted
α-quantiles, q2 (α) = zα + 4b(zα2 ¡ 1)/(3σn3 ) and q3 (α) = (χ2d (α) ¡ d)/(2d)1/2 . Besides, weevaluated empirical quantiles of Ten to assess the accuracy of the proposed approximations
q2 (α) and q3 (α), and the sizes of the proposed tests for some selected parameters.
References[1] Buckley, M.J., Eagleson, G.K., 1988. An approximation to the distribution of quadratic
forms in normal random variables. Austral. J. Statist., 30, 150–159.
[2] Chen, S.X., Qin, Y.L., 2010. A two-sample test for high dimensional data with applications
to gene-set testing. Ann. Statist., 38, 808–835.
[3] Zhang, J.T., 2005. Approximate and asymptotic distributions of chi-squared-type mixtureswith applications. J. Amer. Statist. Assoc., 100, 273–285.

Financial Risk Management with High-Frequency Data
Hiroumi Misaki 1

1

Introduction

Methods for estimating covariance and correlation between multiple asset prices have been
investigated intensively in the field of financial econometrics. In recent decades, daily
covariance and correlations have been estimated directly using high-frequency financial
data or tick-by-tick data. The estimation object is integrated covariance, which is a natural
measure of the covariation of multivariate high-frequency asset prices. Kunitomo andSato [1, 2, 3] proposed a statistical method called the separating information maximum
likelihood (SIML) estimator.
The methods used to estimate the integrated volatility and covariance (subsequently,
correlation) from high-frequency data are called realized measures. In this presentation,
we first compare the accuracy of realised measures using a number of computer simulationsin univariate case. Second, we introduce a brief result when we combine high-frequency
covariance estimation with optimal portfolio selection methods.

2

Robustness of the SIML estimation

We compare the accuracy of realised measures using a number of computer simulations.
We consider a simple realised volatility (RV), a 5-minute RV, a subsampled 5-minute RV,
a two-scale estimator (TS), a realised kernel (RK), a pre-averaging estimator (PA) and aseparating information maximum likelihood estimator (SIML).
We use several non-linear transformation models to obtain the form of the market
microstructure noise. We perform 132 cases of simulation in total.
Our findings are as follows. RV and RV5 are dominated by RV5ss in all cases. If thissimplest class of realised measures is used, subsampling is recommended. RV5ss is heavily
biased when the market microstructure noise or/and round-oﬀ is large. TS, RK and PA
are reasonable in most cases, but we have observed some drawbacks in the case of a large
round-oﬀ with a small noise and a small adjustment with a large noise case. TS is also
problematic in the SSAR model with large noise.
SIML is not biased irrationally in any case; therefore we can conclude that SIML issuﬃciently robust to the form of the market microstructure noise. We have also foundthat SIML is the only realised measure to maintain the consistency in all simulations.
Then, we can conclude that SIML is a reasonable choice in practice when we do not know
the exact form of the market microstructure noise, particularly for assets that are actively
traded.

1

Faculty of Engineering, Information and Systems, University of Tsukuba, Tennodai 1-1-1, Tsukuba
City, Ibaraki 305-8577, JAPAN, hmisaki@risk.tsukuba.ac.jp

3

Portfolio selection

Meucci [4] derived the principal portfolios using principal component analysis. The diversification index is defined as the dispersion of the diversification distribution, based on
its entropy. We maximize the diversification index to obtain the optimal weight w∗ toconstruct the risk diversified portfolio.
We compare the following six methods to construct the risk diversified portfolios:(i) Naive. An equally weighted portfolio, i.e., wi = 1/N .
(ii) OC-125. An analogy of the conventional method of the intraday case; Σ̃t are
obtained by 125-day moving average of outer product of open to close returns.
(iii) SIML-1. We use the previous day’s SIML estimate for asset allocation in thecurrent day: Σ̃t = Σ̂t−1 .
(iv) SIML-5. We use the averages of previous estimates: Σ̃t = (Σ̂t−1 + · · · + Σ̂t−5 )/5
to moderate the estimation errors using a single day.
(v) SIML-22. One-month version of SIML-5: Σ̃t = (Σ̂t−1 + · · · + Σ̂t−22 )/22.
(vi) SIML-0. We use the estimates of the covariance matrix of the same day to carryout portfolio optimization: Σ̃t = Σ̂t .
We selected 10 assets from the stock market. We find that the risk diversified portfolio is feasible when combining with high-frequency covariance estimation. However, theperformances of the portfolios with SIML estimation are not always better than the naive
portfolio. Therefore, a prudent choice of prediction method of covariance matrices isrecommended.

References[1] Kunitomo, N., Sato S.: Separating information maximum likelihood estimation ofrealized volatility and covariance with micro-market noise. Discussion paper CIRJEF-581, Graduate School of Economics, University of Tokyo (2008)[2] Kunitomo, N., Sato S.: The SIML estimation of the integrated volatility of Nikkei-225
futures and hedging coeﬃcients with micro-market noise. Math. Comput. Simulat 8,
1272-1289 (2011)[3] Kunitomo, N., Sato S.: Separating information maximum likelihood estimation of theintegrated volatility and covariance with micro-market noise. N. Am. J. Econ. Financ.
26, 282-309 (2013)[4] Meucci, A.: Managing diversification. Risk, 22(5), 74-7 (2009)














                 

     

          !               ! #     $     % % & &'                                 "            "    !

        % &          "   !

  (     )    ! %,---&      !     

! %,---&        ¾ %1 ,&! )   +  +%,--.&  )    !  *                %)& % & /2   %,-34&

"   

 

%

 





%



&



  











     %    ! 5  / 





 %¡&   7  !

 

  

&










 




  


%




1

%
&

&   

 



¾ %1 &



%1 &





6

,

 

6

,



,
6



&

,         %6 6 &6  ! %,---&      ! 8

%



&

 0







  

  

 



   
%

 %&  % & %  &&

5$   7   %&            







 
 






'   $  

6 &



              ,  

  !

   
 $

 

6  %6

µ

 

 


,

 %&



 


,






 



 
 
 



%

&



  
















 


,





% & 

 



 % & 9










 % &

 


 
    % & )    ! !     ! ! %,-34&!       $           10$ :66;:<.!

$ '!  2   $ =!  2 $ ! %,---&!              
? $ 2! 

5 @@$ )

   A
5 $

7       4$ >.,;>4.!

(! %611,&!                 7$ 4.<;4-.!

! %61,<&! *     B           #0        16$ ,.<;,-6!

6 

High-dimensional analysis of the piecewise deterministic
Markov process for Bayesian inference
Kengo Kamatani (Osaka University, JST CREST)
In this talk, we review the recent results of piecewise deterministic Markov processes forBayesian inference in a high-dimension context. Suppose we wish to sample from a probability
distribution
Πpdxq “ expp´Hpxqqdx px P Rd q

where H : Rd Ñ R is a continuously diﬀerentiable function. For the Bayesian context, thisprobability distribution is the posterior distribution of interest. If we have an i.i.d. sample fromΠ, we can approximate Π-integral of any function f pxq by the law of large numbers
żM
1 ÿf pXm q ÝÑf pxqΠpdxq.
(0.1)
M m“1
Rd

In most of the cases, direct i.i.d. sampling is impossible or computationally very expensive. For
these cases, the Markov chain Monte Carlo method is useful which originated with the classic
paper by Metropolis et al. (1953) almost 70 years ago. The Markov chain Monte Carlo method
is designed to construct an ergodic Markov kernel P which is Π-invariant. If a Markov chain
X1 , X2 , . . . is generaetd from the Markov kernel P then the law of large numbers (0.1) is satisfied.
The Markov chain Monte Carlo is now a gold standard for Bayesian inference.
Recently, its continuous process version, the Markov process Monte Carlo method is ofsubstantial interest for Monte Carlo analysis. Known Markov process Monte Carlo methods relyon an auxiliary variable trick which uses an auxiliary variable v with a probability density ν
on Ξ and considers the joint probability distribution µ :“ Πpdxq b νpdvq as an extended target
distribution on Z “ Rd ˆ Ξ. The original target distribution is a marginal distribution of theextended target distribution. Since Brownian motion does not have an absolutely continuous
path, we can not simulate processes driven by Brownian motion exactly. For our Monte Carloanalysis, exact sampling is necessary. Therefore, the Markov processes of interest should nothave a Brownian part. Known processes consist of a deterministic part and a pure jump part.
These processes are known as the piecewise deterministic Markov processes.
Here we follow Azaïs et al. (2014) for the expression of the piecewise deterministic Markov
processes. The processes are constructed by characteristics pφ, λk , Qk : l “ 1, . . . , Kq. The flow
φ : Z ˆ R Ñ Z is a one-parameter group of homeomorphism, that is, φ is continuous, φp¨, tq isa homeomorphism for each t P R and φpφp¨, sq, tq “ φp¨, s ` tq. For each k “ 1, . . . , K, the jumprate λk : Z Ñ R` determines the jump time of pure jump processes, and Qk is a Markov kernel
şton Z. Let Λk pz, tq “ 0 λk pφpz, sqqds.
The Markov process is defined by the following way. Suppose zp0q “ pxp0q, tp0qq P Z. Let
T1 , . . . , TK be independent processes with PpTk ě tq “ expp´Λk pz, tqq. Let T˚ “ mink“1,...,K Tk .
If Tk “ T˚ , then Z is generated from Qk pφpz, T˚ q, ¨q and set"
φpzp0q, tq for t ă T˚
Xptq “Z
for t “ T˚ .
1

After T˚ , the process evolves in the same way with starting value Z. There are several choicesof characteristics. We introduce three piecewise deterministic Markov processes.
Two popular piecewise deterministic Markov processes use the same flow φ defined by x1 ptq “
vptq and v 1 ptq “ 0. The Zig-Zag sampler proposed by Bierkens et al. (2019) uses d Markov
kernels Q1 , . . . , Qd with d jump rates λ1 , . . . , λd . For each i “ 1, . . . , d, the Markov kernel is adeterministic kernel Qi defined by a map px, vq ÞÑ px, Fi pvqq where Fi is an operator which flips
the sign of the i-th coordinate of x. The jump rate is defined by λi ppx, vqq “ maxt0, Bi Hpxqvi u.
The bouncy particle sampler proposed by Peters and de With (2012), Bouchard-Côtéet al. (2018) uses two Markov kernels Qbounce and Qref with corresponding jump rates λbounceand λref . The kernel Qbounce is a deterministic kernel defined by a map px, vq ÞÑ px, κpx, vqq
wherex∇Hpxq, vyκpx, vq “ v ´ 2∇Hpxq}∇Hpxq}2and λbounce px, vq “ maxt0, x∇Hpxq, vyu. The jump rate λref is a positive constant, and Qref is a µinvariant Markov kernel. For our analysis, for simplicity, we assume Qref ppx, vq, dpy, wqq “ νpdwq.
We review some recent results on asymptotic properties of the above deterministic Markov
processes such as Bierkens et al. (2018), Andrieu et al. (2018), Deligiannidis et al. (2018).

References
Christophe Andrieu, Alain Durmus, Nikolas Nüsken, and Julien Roussel. Hypocoercivity of
Piecewise Deterministic Markov Process-Monte Carlo. arXiv e-prints, art. arXiv:1808.08592,
Aug 2018.
Romain Azaïs, Jean-Baptiste Bardet, Alexandre Génadot, Nathalie Krell, and Pierre-André Zitt.
Piecewise deterministic Markov process—recent results. In Journées MAS 2012, volume 44 of
ESAIM Proc., pages 276–290. EDP Sci., Les Ulis, 2014. doi: 10.1051/proc/201444017.
Joris Bierkens, Kengo Kamatani, and Gareth O. Roberts. High-dimensional scaling limits ofpiecewise deterministic sampling algorithms. arXiv e-prints, art. arXiv:1807.11358, Jul 2018.
Joris Bierkens, Paul Fearnhead, and Gareth Roberts. The zig-zag process and super-eﬃcient
sampling for Bayesian analysis of big data. Ann. Statist., 47(3):1288–1320, 2019. ISSN 00905364. doi: 10.1214/18-AOS1715.
Alexandre Bouchard-Côté, Sebastian J. Vollmer, and Arnaud Doucet. The bouncy particle sampler: A nonreversible rejection-free markov chain monte carlo method. Journal of the American
Statistical Association, 113(522):855–867, 2018. doi: 10.1080/01621459.2017.1294075. URL
https://doi.org/10.1080/01621459.2017.1294075.
George Deligiannidis, Daniel Paulin, Alexandre Bouchard-Côté, and Arnaud Doucet. Randomized Hamiltonian Monte Carlo as Scaling Limit of the Bouncy Particle Sampler and DimensionFree Convergence Rates. arXiv e-prints, art. arXiv:1808.04299, Aug 2018.
N. Metropolis, W. A. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller. Equations ofstate calculations by fast computing machines. Journal of Chemical Physics, 21:1087–1092,
1953.
Elias AJF Peters and G. de With. Rejection-free monte carlo sampling for general potentials.
Physical Review E, 85(2):026703, 2012.

2

A Cauchy family derived by the Möbius transformationsof the sphere
Shogo Kato∗,a and Peter McCullagh b
a The Institute of Statistical Mathematics, Japanb Department of Statistics, University of Chicago, USA

1

A Cauchy family on the sphereThis paper discusses a family of distributions on the unit sphere S d ⊂ Rd+1 with probability
density function()d1 − ρ2Γ{(d + 1)/2}
f (y; µ, ρ) =, y ∈ Sd,(1)
1 + ρ2 − 2ρµT y
2π (d+1)/2
with respect to surface area, where µ ∈ S d is the location parameter, ρ ∈ [0, 1) is the concentration parameter, and S d = {x ∈ Rd+1 ; ∥x∥ = 1} denotes the unit sphere in Rd+1 . The circularcase (d = 1) is well-known as the wrapped Cauchy or circular Cauchy family; see, e.g., Kent
& Tyler (1988) and McCullagh (1996). In this paper, the distribution (1) is called the Cauchy
distribution on the sphere or the spherical Cauchy distribution.
McCullagh (1996) showed that the wrapped Cauchy family is closed under conformal maps
preserving the unit circle which are called the Möbius transformations on the unit circle, andthat there are similar induced transformations on the parameter space. Related results about
the Cauchy family on the real line and on the Euclidean space have been given by McCullagh(1992) and Letac (1986), respectively. To our knowledge, however, there has been no literature
about the association between the Möbius transformations and the spherical Cauchy family(1). Since there have been various statistical applications of the wrapped Cauchy family and/or
the Möbius transformations in directional statistics (McCullagh, 1996; Downs & Mardia, 2002;
Downs, 2003; Jones, 2004; Kato, 2010; Kato & Jones, 2010; Kato & Pewsey, 2015; Uesu el al.,2015), it is potentially useful to consider the Cauchy family on the sphere and its relationshipwith the Möbius transformations.

2

Some properties of a Cauchy family on the sphereThis paper presents some properties of the Cauchy family on the sphere, especially, those related
to the Möbius transformations. The Möbius transformation is defined by}{
1 − ∥ψ∥2(x̃ + ψ) +, x ∈ Rd+1 \ {0, −ψ/∥ψ∥2 }.
(2)MR,ψ (x) = R∥x̃ + ψ∥2where x̃ = x/∥x∥2 , ψ ∈ Rd+1 \S d , and R is a (d + 1) × (d + 1) rotation matrix. Also, we defineMR,ψ (0) = Rψ, MR,ψ (−ψ/∥ψ∥2 ) = ∞ and MR,ψ (∞) = Rψ/∥ψ∥2 . The Möbius transformation (2) is a bijective conformal map which maps R∗

d+1

onto itself, where R

d+1

denotes theAddress for correspondence: Shogo Kato, The Institute of Statistical Mathematics, 10-3 Midori-cho,Tachikawa, Tokyo 190-8562, Japan. E-mail: skato@ism.ac.jp

1

d+1

(d + 1)-dimensional compactified Euclidean space Rd+1 ∪ {∞}. For any ψ ∈ R\ S d , thedtransformation (2) maps the unit sphere S onto itself.
The spherical Cauchy family is closed under the Möbius transformations (2) on the sphere,d+1
and the transformed parameter is given by the Möbius transformation (2) on R . The statistical benefits of this property include: (i) an eﬃcient algorithm for random variate generation;(ii) a simple pivotal statistic for parametric inference; (iii) straightforward calculation of probabilities of a region; (iv) closed form expression for the maximum likelihood estimator for n ≤ 3;
and (v) straightforward calculation of the Fisher information matrix. A method of moments
estimator can be expressed in simple form. A simple algorithm for maximum likelihood estimation is available. The likelihood for the spherical Cauchy is equivalent to that for the t-familywith a certain degree of freedom which is related to the spherical Cauchy via stereographic projection. An asymptotically eﬃcient estimator is presented which our simulation study suggests
outperforms the method of moments estimator and the maximum likelihood estimator in certain
settings. Comparing the densities of the spherical Cauchy and von Mises–Fisher, the spherical
Cauchy density takes greater values around the mode and antimode and smaller values in theother area of the sphere. (See, e.g., Section 9.3.2 of Mardia & Jupp, 1999, for the definition andproperties of the von Mises–Fisher family.) The advantages of the spherical Cauchy over the
von Mises–Fisher in terms of properties include the closure under the Möbius transformations
and the related properties, while the von Mises–Fisher compares favourably with the spherical Cauchy in terms of its membership in the exponential family, straightforward maximum
likelihood estimation and well-developed theory of hypothesis testing.
The preprint version (Kato & McCullagh, 2018) of this paper includes details of the properties of the spherical Cauchy family discussed above.

Acknowledgment
The work of the first author was supported by JSPS KAKENHI Grant Number 17K05379.

References
Downs, T. D. (2003). Spherical regression. Biometrika, 90, 655–668.
Downs, T. D. & Mardia, K. V. (2002). Circular regression. Biometrika, 89, 683–697.
Jones, M. C. (2004). The Möbius distribution on the disc. Ann. Inst. Statist. Math., 56, 733–742.
Kato, S. (2010). A Markov process for circular data. J. R. Statist. Soc. B 72, 655–672.
Kato, S. & Jones, M. C. (2010). A family of distributions on the circle with lines to, and applications
arising from, Möbius transformation. J. Am. Statist. Assoc., 105, 249–262.
Kato, S. & McCullagh, P. (2018). Möbius transformation and a Cauchy family on the sphere. arXiv
preprint arXiv:1510.07679v2 [math.ST].
Kato, S. & Pewsey, A. (2015). A Möbius transformation-induced distribution on the torus.
Biometrika, 102, 359–370.
Kent, J. T. & Tyler, D. E. (1988). Maximum likelihood estimation for the wrapped Cauchy distribution. J. Appl. Statist., 15, 247–254.
Letac, G. (1986). Seul le groupe des similitudes-inversions préserve le type de la loi de Cauchy-conformede Rn pour n > 1. J. Funct. Anal., 68, 43–54.
Mardia, K. V. & Jupp, P. E. (1999). Directional Statistics. Chichester: Wiley.
McCullagh, P. (1992). Conditional inference and Cauchy models. Biometrika, 79, 247–259.
McCullagh, P. (1996). Möbius transformation and Cauchy parameter estimation. Ann. Statist., 24,
787–808.
Uesu, K., Shimizu, K. & SenGupta, A. (2015). A possibly asymmetric multivariate generalizationof the Möbius distribution for directional data. J. Multi. Anal., 134, 146–162.

2

Data Beyond the Euclidean SpaceJörn Schulz
Department of Electrical Engineering and Computer Science,
University of Stavanger, Stavanger, Norway

Complex data such as non-Euclidean or a mixture of Euclidean and nonEuclidean data has gained growing attention recently. However, only few
methods are available to do sensitive statistical inferences on these types ofdata and only little is known about their asymptotic properties. In the following, we assume that the non-Euclidean data lives on a smooth manifold
and in particular we will focus on directional data, i.e. data on the hypersphere Sd “ tx P Rd`1 : xT x “ 1u and data on polyspheres pS2 qd . Examplesof these data types are i.) shape representations including directions such as
skeletal representations that live on Sd1 ˆ pS2 qd2 (Hong et al. (2016); Pizer
et al. (2013); Schulz et al. (2016)), ii.) dihedral angles of protein structures
on pS1 ˆ S1 qd (Eltzner et al. (2017)) or iii.) to analyze temporal sequences ofmolecules on Sd (Dryden et al. (2019)). Especially, in examples i.) and ii.)
we have usually a high dimension low sample size setting, i.e. d " n where
n is the sample size and d is the dimension.
A crucial step in the analysis in all these applications is principal nested
spheres (PNS) (Jung et al. (2012)), a method for decomposition and dimension reduction of directional data on Sd . In opposite to principal component
analysis, PNS is a backward dimension reduction method. In each step,a submanifold of successively lower dimension, containing the largest total
variance, is fitted to the data. A submanifold can be either a small-sphereor a great sphere, i.e. a sphere with radius r ă π{2 or r “ π{2. The choiceof a small or a great sphere is a critical question in the PNS procedure. The
fitting of a small sphere to the data might result in an overfitting, e.g. if
the data is concentrated around a point at Sd . We will discuss a new testingprocedure that outperforms alternative testing methods during a simulation
study and the analysis of skeletal 3D models of hippocampi. The proposed
method is based on a measure of multivariate kurtosis for directional data.
Given a suitable decomposition of the data, statistical inference by hypothesis
testing (Schulz et al. (2016)), classification (Hong et al. (2016)) or clustering(Dryden et al. (2019)) might be performed.
In addition, we will briefly review and discuss some recent works on
asymptotic results within this framework.

1References
Dryden, I. L., Kim, K.-R., Laughton, C. A., and Le, H. (2019), “Principal nested shapespace analysis of molecular dynamics data,” arXiv preprint arXiv:1903.09445.
Eltzner, B., Huckemann, S., and Mardia, K. V. (2017), “Torus Principal Component
Analysis with an Application to RNA Structures,” Annals of Applied Statistics, ISSN
1932-6157 (In Press).
Hong, J., Vicory, J., Schulz, J., Styner, M., Marron, J. S., and Pizer, S. M. (2016), “NonEuclidean classification of medically imaged objects via s-reps,” Medical image analysis,
18, 37–45.
Jung, S., Dryden, I. L., and Marron, J. S. (2012), “Analysis of Principal Nested Spheres,”
Biometrika, 99, 551–568.
Pizer, S. M., Jung, S., Goswami, D., Zhao, X., Chaudhuri, R., Damon, J. N., Huckemann,S., and Marron, J. S. (2013), “Nested Sphere Statistics of Skeletal Models,” in Innovations for Shape Analysis: Models and Algorithms, eds. Breus, M., Bruckstein, A., andMaragos, P., New York: Springer, pp. 93–115.
Schulz, J., Pizer, S. M., Marron, J. S., and Godtliebsen, F. (2016), “Non-linear Hypothesis
Testing of Geometric Object Properties of Shapes Applied to Hippocampi,” Journal of
Mathematical Imaging and Vision, 54, 15–34, issue 1.

2

Ping-Shou Zhong
University of Illinois at Chicago
Title: Change point detection and identification for high dimensional
dependent data
Abstract: High-dimensional functional data appear in practice whena dense number of repeated measurements are taken on a large
number of variables for a relatively small number of experimental
units. The spatial-temporal dependence and high-dimensional natureof the data structure make statistical analysis and computation achallenge. This talk will introduce computationally efficient procedures
to detect and identify change points among covariance matrices fromhigh-dimensional functional data. The change point detection
procedure is presented in the form of a hypothesis test, and theasymptotic distributions of the proposed test statistics are establishedunder an asymptotic framework with “large p, large T and small n”,
where p is data dimension, T is the number of repeatedmeasurements and n is the sample size. We also propose changepoint estimators for both single and multiple change points. Theseestimators are proven to be consistent under a mild set of conditions.
The rate of convergence of the estimator depends on the data
dimension, sample size, number of repeated measurements, andsignal-to-noise ratio. Computation efficiency is carefully studied toaddress the challenges due to the large number of repeatedmeasurements
andhigh-dimensionality.
Simulation
results
demonstrate that the size of the detection procedure is well controlled
at the nominal level, and the locations of multiple change points canaccurately be identified. We apply the proposed approach to findevent boundaries in a continuous movie by identifying change points
among functional connectivity using functional MRI data.

Towards a sparse, scalable, and stably positive definite(inverse) covariance estimator
Joong-Ho (Johann) Won
Department of Statistics, Seoul National University

Abstract: High-dimensional covariance estimation and graphicalmodel selection is a contemporary topic in statistics and machinelearning, and has widespread applications. The problem isnotoriously difficult in high dimensions as the traditional
estimate is not even positive definite, let alone sufficiently
stable. An important line of research is to shrink the spectrum
to yield stable well-conditioned estimators. A separate line ofresearch has considered sparse estimation using nonsmooth
regularization methods and provides interpretable models withfewer parameters. Though an estimator which is both stable andsparse is often desirable in numerous downstream applications,obtaining such estimators is inherently challenging in modern
high-dimensional regimes due to the very different nature of thetwo approaches. In this talk we propose a unifying and scalable
framework which addresses this problem. Our general methodology
takes an arbitrary covariance loss functions (such as the ones
which have been proposed in the literature) and yields estimatesthat are both spectrally regularized and sparse. The frameworkleads to an enriched class of estimators which arecomputationally tractable and enjoy good asymptotic properties.
In addition, when the covariance loss function is orthogonallyinvariant, we further demonstrate that a solution path algorithm
can be derived, involving a series of ordinary differential
equations. The path algorithm is attractive because it provides
the entire family of estimates for all possible values of theregularization parameter, at the same computational cost of asingle estimate with a fixed parameter. An important finding isthat an iterative path algorithm can be devised even when theloss function is not orthogonally invariant, utilizing modern
operator splitting techniques. We illustrate the efficacy of ourapproach on both real and simulated data.
This is a collaboration with Sang-Yun Oh (UC Santa Barbara) andBala Rajaratnam (UC Davis).

A Two-Stage Dimension Reduction Method and
its Applications on Highly Contaminated Image Sets
I-Ping Tu
Institute of Statistical Science, Academia Sinica, Taiwan
Abstract
Principal component analysis (PCA) is arguably the most popular dimension
reduction method for vector type data. When applied on image data, PCA
demands the images to be portrayed as vectors. The resulting computation isheavy because it would solve an eigenvalue problem of a covariance matrix
whose size equals the square of the pixel number. To mitigate the computation
burden, multi-linear PCA that uses column and row basis with a Kronecker
product to compose the matrix structure was proposed, for which the success
was demonstrated on face image sets. However, when we apply MPCA on theparticle images of the single particle cryo-electron microscopy (cryo-EM)
experiments, the results are not satisfying. Here, we propose a dimension
reduction method called Two Stage Dimensional Reduction (2SDR) where wefirst apply MPCA to extract its projection scores, and then apply PCA on thesescores to further reduce the dimension. Tests using single particle cryo-EM
benchmark experimental data sets demonstrate that 2SDR reduce huge
computation costs compared to PCA, and show 2SDR can reconstruct better
quality images than MPCA. Further application of 2SDR on a cryo-EM
micrograph data set significantly reduces the noise to clearly reveal theindividual particles. Remarkably, the de-noised particles boxed out from themicrograph allow subsequent structural analysis to reach a high-quality 3D
density map. This is a joint work with Szu-Chi Chung, Po-Yao Niu, Su-Yun
Huang and Wei-Hau Chang.

Sample covariance matrices from "bad populations"
Jeff Yao
Department of Statistics and Actuarial Science,
The University of Hong Kong

Recent spectral analysis of large covariance matrices is largely basedon the celebrated Marcenko-Pastur law and subsequent applications of
the theory to high-dimensional inference involve central limit theorems
for the corresponding eigenvalue statistics. However it has recently
appeared that there are some important multivariate populations withstrongly dependent coordinates for which the exiting theories do not
apply. High-dimensional mixtures are one of such "bad populations". In
this talk, I will describe this phenomenon and then present somealternative results for the case of high-dimensional mixtures.

Direct estimation of conditional averaging treatmenteﬀect in high dimensions
Shota Katayama
Keio University, JapanThe estimation of conditional average treatment eﬀect (CATE) is a general and fundamental problem in observational studies. Such estimation problem is essential for policy
evaluation, personalized medicine, oﬄine or online marketing and advertising. Usually, toidentify CATE, one requires the strong ignorability condition which says that outcomes andtreatment assignment is independent conditional on covariates. In other words, only thecovariates we collect aﬀect both of outcomes and treatment assignment. If we fail to collectsuch a covariate, the strong ignorability does not hold. Clearly, a large number of covariates
tends to meet the strong ignorability, although it is uncheckable condition from observations. With advances of information technology and database system, it would be plausible
to consider the high dimensional covariates.
In this talk, we consider the estimation of CATE in high dimensions. Following the
Neyman–Rubin’s potential outcome framework (Rubin, 1974, Neyman et al., 1990), assume
that there is a potential outcomes (Yi (0), Yi (1)) for each sample i ∈ {1, 2, . . . , n}. Let Ti ∈{0, 1} be the assignment indicator. Then, Yi (0) ∈ R is the potential outcome when thesample i is assigned to the control (Ti = 0) and Yi (1) ∈ R is the potential outcome when it
is assigned to the treatment (Ti = 1). Assume that we have n independent and identically
distributed examples {(Xi , Ti , Yi (Ti ))}ni=1 where Xi ∈ Rp is the covariates with possibly highdimensions, that is, p ≫ n. Our goal is to estimate the conditional average treatment eﬀect(CATE) given by{}τ ∗ (x) = E Yi (1) − Yi (0) |Xi = x .
To identify the CATE, we assume the following strong ignorability condition.
Assumption 1. {Yi (0), Yi (1)} ⊥⊥ Ti |Xi
Moreover, we assume the linearity for the potential outcomes.
{}{}
Assumption 2. E Yi (0) |Xi = xi = xTi β0∗ and E Yi (1) |Xi = xi = xTi β1∗ .

From Assumption 2, we have τ ∗ (x) = xT (β1∗ − β0∗ ) and we can estimate β1∗ from thetreated examples and can estimate β0∗ from the control examples under Assumption 1. Since
the covariates Xi is high dimension, a natural approach would be applying the Lasso proposed by Tibshirani (1996) for each treated and control examples, i.e.,β̂t = argmin
β∈Rp1∑(Yi − XiT β)2 + λt ∥β∥1 ,2 T =tt = 0, 1,iwhere Yi = Yi (Ti ). Thus, we obtain the estimator of CATE as τ̂ (x) = xT (β̂1 − β̂0 ). However,such the procedure estimate β0∗ and β1∗ separately. The treated (control) outcomes are
predicted by treated (control) covariates. Hence, if x is coming from the distribution of
X |T = 1, then xT β̂1 would be accurate but xT β̂0 be not. Moreover, the non-zero elementsof β̂1 and β̂0 usually do not imply zero elements of β̂1 − β̂0 even when the corresponding
elements of β1∗ − β0∗ are zero.
Our goal is to construct a direct estimation procedure for θ ∗ = β1∗ −β0∗ via the well-known
consequence of the strong ignorability condition, given by]}[ {1 − Ti
Ti−
Xi = x ,τ (x) = E Yie(x) 1 − e(x)∗where e(x) = P(Ti = 1 |Xi = x) is the propensity score function at x. Thus, θ ∗ canbe estimated by regressing the appropriately weighted outcomes on the covariates. The
propensity score is unknown in most cases. An approach to estimate it in high dimensions
may be generalized linear regression with sparse regularization (see, e.g., Fan and Li (2001)
and Van de Geer (2008)), but it may lead to an biased estimator for θ ∗ when the propensity
score function is misspecified.
In this talk, inspired by Athey et al. (2018), a two-step estimation procedure of θ ∗ is
proposed. The first step obtains weightings for outcomes without specifying the propensity
score and then Lasso is applied to the weighted outcomes. The weights are computed by thealternating direction method of multipliers (ADMM) with smoothing technique by Nesterov(2005).

Statistical modeling for electricity load forecasting
Kei Hirose 1,3 and Hiroki Masuda 2
1

Institute of Mathematics for Industry, Kyushu University, 744 Motooka, Nishi-ku, Fukuoka 819-0395,

Japan
2

Faculty of Mathematics, Kyushu University, 744 Motooka, Nishi-ku, Fukuoka 819-0395, Japan

3

RIKEN Center for Advanced Intelligence Project, 1-4-1 Nihonbashi, Chuo-ku, Tokyo 103-0027, Japan

E-mail: hirose@imi.kyushu-u.ac.jp, hiroki@math.kyushu-u.ac.jp

Short and medium-term load forecasting with high accuracy is essential to decision
making in a trade on electricity markets. Among various electricity markets, the dayahead (or spot) market is popular in many electricity exchanges, including Japan ElectricPower Exchange (JEPX) (http://www.jepx.org/english/index.html). In the dayahead market, the contracts for the delivery of electricity on the following day are made,and the transaction is typically carried out in 30 minutes interval; the suppliers must
forecast the loads in 30 minutes interval. In this paper, we consider the problem offorecasting electricity consumption, which can be used for the day-ahead market.
The electricity demand is mainly determined by past electricity consumption andexternal effect, such as weather information. To incorporate the weather information into
the regression model, we may use the weather forecast in 30 minutes interval. In general,
however, the weather forecast information are obtained per not 30 minutes but daily;
for example, only maximum temperature or average humidity is available. Therefore,
ordinary linear regression may not perform well.
In this paper, we introduce statistical modeling which elaborately captures the nonlinear structure of the weather factor with limited weather forecast information. We use thevarying coefficient model with basis expansion, in which the coefficients are assumed tobe different depending on time intervals. The effect of weather variables is expressed as anonlinear function using basis expansions to capture the yearly seasonal effect of weather
information. For interpretation of the estimated model, we eliminate the effect of weather
from the past consumption data. Our regression model turns out to be a linear regression
model, so that we can apply the least-squares estimation.
1

In least-squares estimation, however, we found that the regression coefficients concerning the effect of past electricity consumption often become negative, which makes
the interpretation of the estimated model difficult. To handle this issue, we employ
the non-negative least squares (NNLS; [1], for example) estimation, where we impose aconstraint that the regression coefficients are nonnegative. Furthermore, we employ the
post-selection inference to construct the prediction interval after the model selection basedon [2] and [3].

References[1] D. Chen and R. J. Plemmons. Nonnegativity constraints in numerical analysis. In
The birth of numerical analysis, pages 109–139. World Scientific, 2010.
[2] J. D. Lee and J. E. Taylor. Exact Post Model Selection Inference for Marginal Screening. arXiv.1402.5596, 2014.
[3] J. D. Lee, D. L. Sun, Y. Sun, and J. E. Taylor. Exact post-selection inference, withapplication to the lasso. The Annals of Statistics, 44(3):907–927, June 2016.

2

Consistency of the objective general index
in high dimensional settings
Takuma Bando∗, Tomonari Sei∗ and Kazuyoshi Yata†

1

Introduction

Rankings are often determined by multivariate data. For example, the world
university ranking provided by [3] is based on five attributes of universities:
teaching, research, citations, industry income and international outlook. Fora happiness index of prefectures in Japan [2], 65 attributes are used to makea ranking of 47 prefectures. In heptathlon of athletics, the scores of seven
events are unified into an overall score. These rankings are, after sometransformations, based on a weighted sum of variables.
We focus on the weights. In [1], an objective weight is proposed via diagonal scaling of the sample covariance matrix. The resultant index called theobjective general index (OGI) has positive correlation with all the variables
and is invariant with respect to scale transformation of the data.
In some applications like the happiness ranking mentioned above, thenumber of variables is often large and comparable with the sample size. In
other words, we have to deal with high-dimensional data for ranking. If
we use the objective general index for such a data, a reliable estimator will
be required. The aim of this paper is to study consistency of the weight
determined from a random sample.

2

Problem setting and main results

The objective general index is one of possible general indices for multivariate
data. We consider a weighted sum w⊤ x of an observation x ∈ Rp as a∗†

The University of Tokyo
University of Tsukuba

general index, where w ∈ Rp is a weight vector. Each variate xi is assumed
to have a meaning that “larger is better” without loss of generality. Then
it is natural to suppose that every coordinate of w is positive.
Let n > p and consider a random sample x(1) , . . . , x(n) ∈ Rp according
to the multivariate normal distribution N (0, Σ) with covariance matrix Σ ∈∑Rp×p . Denote the sample covariance matrix by S = n−1 nt=1 x(t) x⊤(t) .
Definition 1 ([1]). The objective weight w is defined by a solution of
Σw =

1
,
w

w ∈ Rp>0 .

(1)Similarly, the sample objective weight ŵ is defined by
S ŵ =

1,ŵ

ŵ ∈ Rp>0 .

(2)

The weighted sum w⊤ x of an observation x ∈ Rp using the objective weight
w is called the objective general index (OGI).
We consider a high-dimensional setting in that the dimension p grows
with the sample size n. Denote the entries of Σ as (σij )pi,j=1 .
Theorem 1. Suppose that Σ1 = 1. Then there exists a constant C > 0such that(nCε2P(∥ŵ − 1∥ ≥ ε) ≤ 4p exp −(maxi σii )p2)(3)

for any ε > 0 and any n ≥ n0 with some n0 = n0 (ε). In particular, if
maxi σii = O(1) and (p2 log p)/n = o(1) as n → ∞, then ŵ is strongly
consistent in the sense that ∥ŵ − 1∥2 converges to 0 almost surely.

References[1] Sei, T. (2016). An objective general index for multivariate ordered data, J.
Multivariate Anal., 147, 247–264.
[2] Terashima, J., Japan Research Institute / Nihon Unisys Ltd. (eds.) (2016)
Happiness ranking of all the 47 prefectures in Japan, 2016 (in Japanese), Toyo
Keizai.
[3] Times Higher Education, World University Rankings.
https://www.timeshighereducation.com/