# High-order Uncertain Differential Equation and Its Application to Nuclear
Reactors

Hao Li, Yuqian Wang∗ School of Mathematics, Renmin University of China,
Beijing 100872, China. Email<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

High-order uncertain differential equation (HUDE) was introduced in
literature. But the present method to solve a HUDE is incorrect. In this
paper, we will rigorously prove some comparion theorems of high-order
differential equations, and present a method to solve a family of HUDE,
including parameter estimation and hypothesis test. Then an application to
nuclear reactor kinetics is given to illustrate the method.

Key words: Uncertainty theory; Uncertain differential equation; High-order
uncertain differential equation; Nuclear reactor kinetics.

##  1 Introduction

Uncertainty theory, founded by [5], has been developed into an axiomatic
mathematical theory. Among the many theoretical branches of uncertainty
theory, uncertain statistics has been the most popular and cutting-edge
theoretical branch for the last few years. There are three important methods
in uncertain statistics: uncertain time series analysis, uncertain regression
analysis and uncertain differential equations. They have drawn attention of
researchers with different background. Many achievements and results have been
made in theory and applications.

Uncertain differential equation was first introduced by [6]. For those
uncertain differential equations with no analytical solutions, [17] presented
a useful formula to calculate the inverse uncertainty distribution of the
solution to an uncertain differential equation in terms of $\alpha$-path. This
formula is known as ”Yao-Chen Formula”. [21] discussed the uncertain partial
differential equations. [18] discussed the partial derivatives of uncertain
field, and gave the integral form of an uncertain partial differential
equations.

Uncertain differential equations has been widely applied to many fields, such
as chemical reaction ([13]), pharmacokinetics ([11]), epidemic spread ([4]),
gas price ([12]), China’s population([14]), China’s birth rates([20]).
Recently, some practical analysis in finance verify that, compared with
stochastic differential equations, uncertain differential equations are more
suitable to fit the data, for example, Alibaba stock price ([8]), currency
exchange rate ([19]), interest rate ([15]). When applying uncertain
differential equations to practical problems, there are two core problems: how
to estimate unknown parameters in an uncertain differential equation based on
the observed data, and how to test the fitness of an uncertain differential
equation. [8] introduced the concept of residuals of uncertain differential
equations, and developed the method of moments estimation. Soon, [9] explored
a modified maximum likelihood estimation, and [10] presented the least squares
estimation. For evaluating an uncertain differential equation’s goodness of
fit, [19] intorduced the hypothesis test.

For complex dynamic systems, such as spring vibration, pendulum swing, RLC
circuit, nuclear reactor kinetics, high-order uncertain differential equations
are required to characterize them. [16] initially proposed the HUDE. However,
in the process of research, we found that there is a mistake in [16] that the
theorem of $\alpha$-path of HUDE is not strictly proved. Therefore, the mathod
of solving a high-order uncertian differentail equation in [16] is wrong.

In this paper, we focus on high-order uncertain differential equations. We
will prove some comparion theorems of high-order ordinary differential
equation, and rigorously prove a theorem of $\alpha$-path of HUDE. Then we
present a method to solve a HUDE, including parameter estimation and
hypothesis test. As an application, the nuclear reactor kinetics under
uncertain circumstance is discussed. The remainder of this paper is organized
as follows. In Section 2, we will prove comparison theorems and give the
solution of a HUDE. In Section 3, we will introduce the concept of residual of
a HUDE. We will discuss parameter estimation and hypothesis test in Section 4
and Section 5, respectively. In Section 6, the model of HUDE will be applied
to nuclear reactor kinetics. Finally, some conclusions will be made in Section
7.

##  2 Solution of a High-order Uncertain Differential Equation

An uncertain differential equation

$\frac{\mathrm{d}^{n}X_{t}}{\mathrm{d}t^{n}}=f\bigg{(}t,X_{t},\cdots,\frac{\mathrm{d}^{n-1}X_{t}}{\mathrm{d}t^{n-1}}\bigg{)}+\sum_{i=1}^{m}g_{i}\bigg{(}t,X_{t},\cdots,\frac{\mathrm{d}^{n-1}X_{t}}{\mathrm{d}t^{n-1}}\bigg{)}\frac{\mathrm{d}C_{it}}{\mathrm{d}t}$
(1)

is called a high-order uncertain differential equation, where $f$ and $g_{i}\
(i=1,2,\cdots,m)$ are continuous functions, and $C_{1t},C_{2t},\cdots,C_{mt}$
are independent Liu processes. Since most of the high-order uncertain
differential equations cannot be solved analytically, the uncertainty
distribution of the solution cannot be determined. So it is significantly
important to figure out the inverse uncertainty distribution of the solution.

For a high-order uncertain differential equation (1), let $X_{t}^{\alpha}$ be
the solution of the corresponding ordinary differential equation:

$\frac{\mathrm{d}^{n}X_{t}^{\alpha}}{\mathrm{d}t^{n}}=f\bigg{(}t,X_{t}^{\alpha},\cdots,\frac{\mathrm{d}^{n-1}X_{t}^{\alpha}}{\mathrm{d}t^{n-1}}\bigg{)}+\sum_{i=1}^{m}\bigg{|}g_{i}\bigg{(}t,X_{t}^{\alpha},\cdots,\frac{\mathrm{d}^{n-1}X_{t}^{\alpha}}{\mathrm{d}t^{n-1}}\bigg{)}\bigg{|}\Phi^{-1}(\alpha)$
(2)

where

$\Phi^{-1}(\alpha)=\frac{\sqrt{3}}{\pi}\mathrm{ln}\frac{\alpha}{1-\alpha},\qquad
0<\alpha<1.$

In [16], $X_{t}^{\alpha}$ is proved to be the $\alpha$-path of $X_{t}$, where
$X_{t}$ is the solution of (1). The relevant conclusion is listed below.

Conclusion:([16], Page 144, Theorem 10.1) The solution $X_{t}$ of a high-order
uncertain differential equation

$\frac{\mathrm{d}^{n}X_{t}}{\mathrm{d}t^{n}}=f\bigg{(}t,X_{t},\cdots,\frac{\mathrm{d}^{n-1}X_{t}}{\mathrm{d}t^{n-1}}\bigg{)}+g\bigg{(}t,X_{t},\cdots,\frac{\mathrm{d}^{n-1}X_{t}}{\mathrm{d}t^{n-1}}\bigg{)}\frac{\mathrm{d}C_{t}}{\mathrm{d}t}$

is a contour process with an $\alpha$-path $X_{t}^{\alpha}$ that solves the
corresponding high-order ordinary differential equation

$\frac{\mathrm{d}^{n}X_{t}^{\alpha}}{\mathrm{d}t^{n}}=f\bigg{(}t,X_{t}^{\alpha},\cdots,\frac{\mathrm{d}^{n-1}X_{t}^{\alpha}}{\mathrm{d}t^{n-1}}\bigg{)}+\bigg{|}g\bigg{(}t,X_{t}^{\alpha},\cdots,\frac{\mathrm{d}^{n-1}X_{t}^{\alpha}}{\mathrm{d}t^{n-1}}\bigg{)}\bigg{|}\Phi^{-1}(\alpha)$

where

$\Phi^{-1}(\alpha)=\frac{\sqrt{3}}{\pi}\mathrm{ln}\frac{\alpha}{1-\alpha},\qquad
0<\alpha<1$

is the inverse uncertainty distribution of standard normal uncertain
variables. In other words,

$\mbox{\euler M}\big{\\{}X_{t}\leq X_{t}^{\alpha},\forall t\big{\\}}=\alpha,$
$\mbox{\euler M}\big{\\{}X_{t}\geq X_{t}^{\alpha},\forall
t\big{\\}}=1-\alpha.$

However, there is a drawback in its proof. Here is a counterexample.

######  Example 2.1.

Consider the following uncertain differential equation when $t\geq 0$

$\left\\{\begin{array}[]{ll}\vspace{2mm}&\displaystyle\frac{\mathrm{d}^{2}X_{t}}{\mathrm{d}t^{2}}=-X_{t}+e^{-t}\frac{\mathrm{d}C_{t}}{\mathrm{d}t}\\\
\vspace{2mm}&\displaystyle X_{0}=a\\\
&\displaystyle\frac{\mathrm{d}X_{t}}{\mathrm{d}t}\bigg{|}_{t=0}=b\end{array}\right.\hskip
30.0pt$ (3)

where $a$ and $b$ are constants. Suppose that $\Psi_{t}^{-1}(\alpha)$ is the
inverse uncertainty distribution of $X_{t}$. If the conclusion is correct,
then we have

$\Psi_{t}^{-1}(\alpha)=X_{t}^{\alpha}.$ (4)

Note that $\Psi_{t}^{-1}(\alpha)$ is an increasing function with respect to
$\alpha$. Choose $\alpha_{1},\alpha_{2}$ such that
$0<\alpha_{1}<\alpha_{2}<1$. Let $p_{1}=\Phi^{-1}(\alpha_{1})$,
$p_{2}=\Phi^{-1}(\alpha_{2})$, and then $p_{1}<p_{2}$. Assume
$X_{t}^{\alpha_{i}}(i=1,2)$ is the solution of

$\left\\{\begin{array}[]{ll}\vspace{2mm}&\displaystyle\frac{\mathrm{d}^{2}X_{t}}{\mathrm{d}t^{2}}=-X_{t}+p_{i}e^{-t}\\\
\vspace{2mm}&\displaystyle X_{0}=a\\\
&\displaystyle\frac{\mathrm{d}X_{t}}{\mathrm{d}t}\bigg{|}_{t=0}=b.\end{array}\right.\hskip
30.0pt$ (5)

It is easy to check that

$X_{t}^{\alpha_{i}}=\Big{(}a-\frac{p_{i}}{2}\Big{)}\cos{t}+\Big{(}b+\frac{p_{i}}{2}\Big{)}\sin{t}+\frac{p_{i}}{2}e^{-t},\quad
i=1,2.$

The graphs of $X_{t}^{\alpha}$ of some different values of $\alpha$ are shown
in Fig. 1. Then

$X_{t}^{\alpha_{1}}-X_{t}^{\alpha_{2}}=\frac{1}{2}(p_{1}-p_{2})\bigg{(}\sqrt{2}\sin\Big{(}t-\frac{\pi}{4}\Big{)}+e^{-t}\bigg{)}.$
(6)

Figure 1: The graphs of $X_{t}^{\alpha}$ of some different values of $\alpha$

So

$X_{t}^{\alpha_{1}}-X_{t}^{\alpha_{2}}>0$

for

$t\in\bigg{(}\displaystyle-\frac{3\pi}{4}+2k\pi,\frac{\pi}{4}+2k\pi\bigg{)}\
(k\in\mathbb{N}^{+}).$

Thus, by Eq.(4),

$\Psi_{t}^{-1}(\alpha_{1})=X_{t}^{\alpha_{1}}>X_{t}^{\alpha_{2}}=\Psi_{t}^{-1}(\alpha_{2}).$

This is a contradiction to the fact that $\Psi_{t}^{-1}(\alpha)$ is an
increasing function. Therefore, this conclusion is incorrect.

Next, we will prove two comparison theorems of ordinary differential equation,
and give a sufficient condition for $X_{t}^{\alpha}$ to be the $\alpha$-path
of $X_{t}$.

######  Theorem 2.1.

Let $f\big{(}t,z_{0},z_{1},\cdots,z_{n-1}\big{)}$ and
$g\big{(}t,z_{0},z_{1},\cdots,z_{n-1}\big{)}$ be two functions on
$D=[t_{0},a]\\!\times\\!\mathbb{D}_{0}\subset\mathbb{R}^{n\\!+\\!1}$ (where
$\mathbb{D}_{0}\ \mbox{is an n-dimensional bounded region})$ satisfying local
Lipschitz conditions in $z_{0},z_{1},\cdots,z_{n-1}$. Assume that
$f\big{(}t,z_{0},z_{1},$
$\cdots,z_{n-1}\big{)}<g\big{(}t,z_{0},z_{1},\cdots,z_{n-1}\big{)}$ on $D$,
and that at least one of $f\big{(}t,z_{0},z_{1},\cdots,z_{n-1}\big{)}$ and
$g\big{(}t,z_{0},z_{1},\cdots,z_{n-1}\big{)}$ is a monotonically increasing
function with respect to $z_{0},z_{1},\cdots,$ $z_{n-2}$. If $\psi(t)$ and
$\Psi(t)$ are solutions of

$\left\\{\begin{aligned}
&\frac{\mathrm{d}^{n}X_{t}}{\mathrm{d}t^{n}}=f\bigg{(}t,X_{t},\cdots,\frac{\mathrm{d}^{n\\!-\\!1}X_{t}}{\mathrm{d}t^{n\\!-\\!1}}\bigg{)}\\\
&X(t_{0})=y_{0},\frac{\mathrm{d}X_{t}}{\mathrm{d}t}\bigg{|}_{t=t_{0}}=y^{\prime}_{0},\cdots,\frac{\mathrm{d}^{n-1}X_{t}}{\mathrm{d}t^{n-1}}\bigg{|}_{t=t_{0}}=y_{0}^{(n-1)}\end{aligned}\right.$

and

$\left\\{\begin{aligned}
&\frac{\mathrm{d}^{n}X_{t}}{\mathrm{d}t^{n}}=g\bigg{(}t,X_{t},\cdots,\frac{\mathrm{d}^{n\\!-\\!1}X_{t}}{\mathrm{d}t^{n\\!-\\!1}}\bigg{)}\\\
&X(t_{0})=y_{0},\frac{\mathrm{d}X_{t}}{\mathrm{d}t}\bigg{|}_{t=t_{0}}=y^{\prime}_{0},\cdots,\frac{\mathrm{d}^{n-1}X_{t}}{\mathrm{d}t^{n-1}}\bigg{|}_{t=t_{0}}=y_{0}^{(n-1)},\end{aligned}\right.$

respectively, then

$\psi(t)<\Psi(t),\ \forall t\in[t_{0},a].$

$\mathbf{Proof.}$ Let

$\phi(t)=\Psi(t)-\psi(t),\quad t\in[t_{0},a].$

Then $\phi(t_{0})=\phi^{\prime}(t_{0})=\cdots=\phi^{(n-1)}(t_{0})=0,$ and
$\phi^{(n)}(t_{0})>0$. By the sign-preserving property of derivatives, there
exists $\delta>t_{0}$ such that

$\phi(t)>0,\ \phi^{\prime}(t)>0,\cdots,\phi^{(n-1)}(t)>0,\ \forall
t\in(t_{0},\delta).$ (7)

By way of contradiction, suppose that $\phi(t)$ is not always positive when
$t\in[t_{0},a]$. Let

$t_{1}=\min\big{\\{}t~{}|~{}\phi(t)=0,\ t\in[t_{0},a]\big{\\}}.$

Then

$\delta\leq t_{1},\mbox{ and }\phi(t)>0,\ \forall t\in(t_{0},t_{1}).$ (8)

By the choice of $t_{1}$ and the above inequalities,

$\phi^{\prime}(t_{1})\leq 0.$ (9)

Next, we define $t_{2},t_{3},\cdots,t_{n-1}$ and $t_{n}$ one by one in the
following way. For $i\ (2\leq i\leq n)$, suppose $t_{i-1}$ is defined such
that

$\delta\leq t_{i-1},\quad\phi^{(i-1)}(t_{i-1})\leq 0$

and

$\phi(t)>0,\phi^{\prime}(t)>0,\cdots,\phi^{(i-2)}(t)>0,\ \forall
t\in(t_{0},t_{i-1}).$

Let

$t_{i}=\min\Big{\\{}t~{}\big{|}~{}\phi^{(i-1)}(t)=0,\ t\in[t_{0},a]\Big{\\}}.$
(10)

Then $\delta\leq t_{i}\leq t_{i-1}$ and

$\phi(t)>0,\phi^{\prime}(t)>0,\cdots,\phi^{(i-1)}(t)>0,\ \forall
t\in(t_{0},t_{i}).$ (11)

By the choice of $t_{i}$ and the above inequalities,

$\phi^{(i)}(t_{i})\leq 0.$ (12)

Repeat this process until we have $t_{1},t_{2},\cdots,t_{n}$. By the
definition of $t_{i}$ ($1\leq i\leq n$), the sequence
$\\{t_{1},t_{2},\cdots,t_{n}\\}$ is decreasing with a lower bound $\delta$
(Fig. 2).

Figure 2: $\\{t_{1},t_{2},\cdots,t_{n}\\}$ is decreasing with a lower bound
$\delta$

Let $k$ be the minimum number such that $t_{k}=t_{n}$. Then
$t_{k-1}>t_{k}=t_{k+1}=\cdots=t_{n}$. By the choice of $t_{n}$,

$\phi(t)>0,\phi^{\prime}(t)>0,\cdots,\phi^{(n-1)}(t)>0,\ \forall
t\in(t_{0},t_{k})$

and

$\phi^{(n)}(t_{k})\leq 0.$ (13)

As $t_{k}=t_{k+1}=\dots=t_{n}$, $\phi^{(i)}(t_{k})=0$ for $k-1\leq i\leq n-1$
by Eq.(10). That is,

$\Psi^{(i)}(t_{k})=\psi^{(i)}(t_{k}),\mbox{ for }k-1\leq i\leq n-1.$ (14)

As $t_{k}<t_{k-1}$, by Eq.(11), $\phi^{(i)}(t_{k})>0$ for $0\leq i\leq k-2$,
i.e.

$\Psi^{(i)}(t_{k})>\psi^{(i)}(t_{k}),\mbox{ for }0\leq i\leq k-2.$ (15)

If $g\big{(}t,z_{0},z_{1},\cdots,z_{n-1}\big{)}$ is a monotonically increasing
function with respect to $z_{0},\cdots,$ $z_{n-2}$, by Eq. (14) and (15),

$\displaystyle\phi^{(n)}(t_{k})$
$\displaystyle=\Psi^{(n)}(t_{k})-\psi^{(n)}(t_{k})$
$\displaystyle=g\Big{(}t_{k},\Psi(t_{k}),\cdots,\Psi^{(n-1)}(t_{k})\Big{)}-f\Big{(}t_{k},\psi(t_{k}),\cdots,\psi^{(n-1)}(t_{k})\Big{)}$
$\displaystyle\geq
g\Big{(}t_{k},\psi(t_{k}),\cdots,\psi^{(n-1)}(t_{k})\Big{)}-f\Big{(}t_{k},\psi(t_{k}),\cdots,\psi^{(n-1)}(t_{k})\Big{)}$
$\displaystyle>0.$

If $f\big{(}t,z_{0},z_{1},\cdots,z_{n-1}\big{)}$ is a monotonically increasing
function with respect to $z_{0},\cdots,$ $z_{n-2}$, by Eq. (14) and (15),

$\displaystyle\phi^{(n)}(t_{k})$
$\displaystyle=g\Big{(}t_{k},\Psi(t_{k}),\cdots,\Psi^{(n-1)}(t_{k})\Big{)}-f\Big{(}t_{k},\psi(t_{k}),\cdots,\psi^{(n-1)}(t_{k})\Big{)}$
$\displaystyle\geq
g\Big{(}t_{k},\Psi(t_{k}),\cdots,\Psi^{(n-1)}(t_{k})\Big{)}-f\Big{(}t_{k},\Psi(t_{k}),\cdots,\Psi^{(n-1)}(t_{k})\Big{)}$
$\displaystyle>0.$

This contradicts to Eq.(13). Therefore this assumption is not valid and the
conclusion $\phi(t)>0$ for $t\in[t_{0},a]$ is true, i.e.

$\psi(t)<\Psi(t),\ \forall t\in[t_{0},a].$

The theorem is proved.

######  Theorem 2.2.

Let $f\big{(}t,z_{0},z_{1},\cdots,z_{n-1}\big{)}$ and
$g\big{(}t,z_{0},z_{1},\cdots,z_{n-1}\big{)}$ be two functions on
$D=[t_{0},a]\\!\times\\!\mathbb{D}_{0}\subset\mathbb{R}^{n\\!+\\!1}$ (where
$\mathbb{D}_{0}\ \mbox{is an n-dimensional bounded region})$ satisfying local
Lipschitz conditions in $z_{0},z_{1},\cdots,z_{n-1}$. Assume that
$f\big{(}t,z_{0},z_{1},$ $\cdots,z_{n-1}\big{)}\leq
g\big{(}t,z_{0},z_{1},\cdots,z_{n-1}\big{)}$ on $D$, and that at least one of
$f\big{(}t,z_{0},z_{1},\cdots,z_{n-1}\big{)}$ and
$g\big{(}t,z_{0},z_{1},\cdots,z_{n-1}\big{)}$ is a monotonically increasing
function with respect to $z_{0},z_{1},\cdots,$ $z_{n-2}$. If $\psi(t)$ and
$\Psi(t)$ are solutions of

$\left\\{\begin{aligned}
&\frac{\mathrm{d}^{n}X_{t}}{\mathrm{d}t^{n}}=f\bigg{(}t,X_{t},\cdots,\frac{\mathrm{d}^{n\\!-\\!1}X_{t}}{\mathrm{d}t^{n\\!-\\!1}}\bigg{)}\\\
&X(t_{0})=y_{0},\frac{\mathrm{d}X_{t}}{\mathrm{d}t}\bigg{|}_{t=t_{0}}=y^{\prime}_{0},\cdots,\frac{\mathrm{d}^{n-1}X_{t}}{\mathrm{d}t^{n-1}}\bigg{|}_{t=t_{0}}=y_{0}^{(n-1)}\end{aligned}\right.$

and

$\left\\{\begin{aligned}
&\frac{\mathrm{d}^{n}X_{t}}{\mathrm{d}t^{n}}=g\bigg{(}t,X_{t},\cdots,\frac{\mathrm{d}^{n\\!-\\!1}X_{t}}{\mathrm{d}t^{n\\!-\\!1}}\bigg{)}\\\
&X(t_{0})=y_{0},\frac{\mathrm{d}X_{t}}{\mathrm{d}t}\bigg{|}_{t=t_{0}}=y^{\prime}_{0},\cdots,\frac{\mathrm{d}^{n-1}X_{t}}{\mathrm{d}t^{n-1}}\bigg{|}_{t=t_{0}}=y_{0}^{(n-1)},\end{aligned}\right.$

respectively, then

$\psi(t)\leq\Psi(t),\ \forall t\in[t_{0},a].$

$\mathbf{Proof.}$ Let $\\{\varepsilon_{i}\\}\ (i=1,2,\cdots)$ be a
monotonically decreasing sequence of positive numbers such that
$\lim\limits_{i\rightarrow\infty}\varepsilon_{i}=0$. Due to

$f\big{(}t,z_{0},z_{1},\cdots,z_{n-1}\big{)}\leq
g\big{(}t,z_{0},z_{1},\cdots,z_{n-1}\big{)}$

for $\big{(}t,z_{0},z_{1},\cdots,z_{n-1}\big{)}\in D$, we have

$f\big{(}t,z_{0},z_{1},\cdots,z_{n-1}\big{)}-\varepsilon_{i}<g\big{(}t,z_{0},z_{1},\cdots,z_{n-1}\big{)}$

for $\big{(}t,z_{0},z_{1},\cdots,z_{n-1}\big{)}\in D$, where $i=1,2,\cdots$.

Consider equations

$\left\\{\begin{aligned}
&\frac{\mathrm{d}^{n}X_{t}}{\mathrm{d}t^{n}}=f\bigg{(}t,X_{t},\cdots,\frac{\mathrm{d}^{n\\!-\\!1}X_{t}}{\mathrm{d}t^{n\\!-\\!1}}\bigg{)}-\varepsilon_{i}\\\
&X(t_{0})=y_{0},\frac{\mathrm{d}X_{t}}{\mathrm{d}t}\bigg{|}_{t=t_{0}}=y^{\prime}_{0},\cdots,\frac{\mathrm{d}^{n-1}X_{t}}{\mathrm{d}t^{n-1}}\bigg{|}_{t=t_{0}}=y_{0}^{(n-1)}.\end{aligned}\right.$
(16)

From the existence and uniqueness theorem of solution, the initial value
problem (16) has exactly one solution $\psi_{i}(t)\ (i=1,2,\cdots)$ in the
interval $t_{0}\leq t\leq a$. According to Theorem 2.1, we can get

$\psi_{i}(t)<\Psi(t)\ (i=1,2,\cdots),\ \forall t\in[t_{0},a].$

Note that

$\displaystyle\big{|}\psi_{i}(t_{1})-\psi_{i}(t_{2})\big{|}$
$\displaystyle\leq\int_{t_{1}}^{t_{2}}\idotsint_{\mathbb{D}_{0}}\bigg{|}f\bigg{(}t,X_{t},\cdots,\frac{\mathrm{d}^{n\\!-\\!1}X_{t}}{\mathrm{d}t^{n\\!-\\!1}}\bigg{)}-\varepsilon_{i}\bigg{|}\,dX_{t}\dots
dt$
$\displaystyle\leq\int_{t_{1}}^{t_{2}}\bigg{(}\idotsint_{\mathbb{D}_{0}}M\,dX_{t}\dots\bigg{)}dt$
$\displaystyle\leq
M\cdot\left\|\mathbb{D}_{0}\right\|\cdot\left|t_{2}-t_{1}\right|$

$\big{(}t_{1},t_{2}\in[t_{0},a],i\\!=\\!1,2,\\!\cdots\\!\big{)}$, where

$M\\!=\\!\max\limits_{\big{(}t,X_{t},\cdots,\frac{\mathrm{d}^{n\\!-\\!1}X_{t}}{\mathrm{d}t^{n\\!-\\!1}}\big{)}\in
D}\left|f\bigg{(}t,X_{t},\\!\cdots\\!,\frac{\mathrm{d}^{n\\!-\\!1}X_{t}}{\mathrm{d}t^{n\\!-\\!1}}\bigg{)}\\!-\\!\varepsilon_{i}\right|.$

Therefore, $\psi_{i}(t)\ (i=1,2,\cdots)$ is uniformly bounded and equally
continuous on $[t_{0},a]$ based on local Lipschitz condition. It follows from
Ascoli Lemma that $\psi_{i}(t)\ (i=1,2,\cdots)$ has uniformly convergent
subsequences on region $D$. Then

$\lim_{i\rightarrow\infty}\psi_{i}(t)=\psi(t).$

Consequently,

$\psi(t)=\lim_{i\rightarrow\infty}\psi_{i}(t)\leq\Psi(t),\ \forall
t\in[t_{0},a].$

The theorem is proved.

Next, we will prove the theorem of $\alpha$-path for HUDE.

######  Theorem 2.3.

Let $X_{t}$ and $X_{t}^{\alpha}$ be the solution of

$\frac{\mathrm{d}^{n}X_{t}}{\mathrm{d}t^{n}}=f\bigg{(}t,X_{t},\frac{\mathrm{d}X_{t}}{\mathrm{d}t},\cdots,\frac{\mathrm{d}^{n-1}X_{t}}{\mathrm{d}t^{n-1}}\bigg{)}+\sum_{i=1}^{m}g_{i}\bigg{(}t,X_{t},\frac{\mathrm{d}X_{t}}{\mathrm{d}t},\cdots,\frac{\mathrm{d}^{n-1}X_{t}}{\mathrm{d}t^{n-1}}\bigg{)}\frac{\mathrm{d}C_{it}}{\mathrm{d}t}$

and

$\frac{\mathrm{d}^{n}X_{t}^{\alpha}}{\mathrm{d}t^{n}}=f\bigg{(}t,X_{t}^{\alpha},\cdots,\frac{\mathrm{d}^{n-1}X_{t}^{\alpha}}{\mathrm{d}t^{n-1}}\bigg{)}+\sum_{i=1}^{m}\bigg{|}g_{i}\bigg{(}t,X_{t}^{\alpha},\cdots,\frac{\mathrm{d}^{n-1}X_{t}^{\alpha}}{\mathrm{d}t^{n-1}}\bigg{)}\bigg{|}\Phi^{-1}(\alpha)$

respectively, where

$\Phi^{-1}(\alpha)=\frac{\sqrt{3}}{\pi}\mathrm{ln}\frac{\alpha}{1-\alpha},\qquad
0<\alpha<1.$

If

$f\bigg{(}t,X_{t}^{\alpha},\frac{\mathrm{d}X_{t}^{\alpha}}{\mathrm{d}t},\cdots,\frac{\mathrm{d}^{n-1}X_{t}^{\alpha}}{\mathrm{d}t^{n-1}}\bigg{)}+\sum\limits_{i=1}\limits^{m}\bigg{|}g_{i}\bigg{(}t,X_{t}^{\alpha},\frac{\mathrm{d}X_{t}^{\alpha}}{\mathrm{d}t},\cdots,\frac{\mathrm{d}^{n-1}X_{t}^{\alpha}}{\mathrm{d}t^{n-1}}\bigg{)}\bigg{|}\Phi^{-1}(\alpha)$

is a monotonically increasing function with respect to

$X_{t}^{\alpha},\frac{\mathrm{d}X_{t}^{\alpha}}{\mathrm{d}t},\cdots,\frac{\mathrm{d}^{n-2}X_{t}^{\alpha}}{\mathrm{d}t^{n-2}},$

then $X_{t}^{\alpha}$ is the $\alpha$-path of $X_{t}$, i.e.,

$\displaystyle\mbox{\euler M}\big{\\{}X_{t}\leq X_{t}^{\alpha},\forall
t\big{\\}}=\alpha,$ $\displaystyle\mbox{\euler
M}\big{\\{}X_{t}>X_{t}^{\alpha},\forall t\big{\\}}=1-\alpha.$

$\mathbf{Proof.}$ Given $\alpha\in(0,1)$, for each $X_{t}^{\alpha}$, we can
divide it into two parts,

$T_{i}^{+}=\bigg{\\{}t\bigg{|}g_{i}\bigg{(}t,X_{t}^{\alpha},\cdots,\frac{\mathrm{d}^{n-1}X_{t}^{\alpha}}{\mathrm{d}t^{n-1}}\bigg{)}\geq
0\bigg{\\}},$
$T_{i}^{-}=\bigg{\\{}t\bigg{|}g_{i}\bigg{(}t,X_{t}^{\alpha},\cdots,\frac{\mathrm{d}^{n-1}X_{t}^{\alpha}}{\mathrm{d}t^{n-1}}\bigg{)}<0\bigg{\\}},$

$i=1,2,\cdots,m$. It is obvious that $T_{i}^{+}\cap T_{i}^{-}=\varnothing$ and
$T_{i}^{+}\cup T_{i}^{-}=[0,+\infty)$ for each $i$, $1\leq i\leq m$.

Next, we define

$\displaystyle\Lambda_{i1}^{+}=\left\\{\gamma\bigg{|}\frac{\mathrm{d}C_{it}(\gamma)}{\mathrm{d}t}\leq\Phi^{-1}(\alpha)\
\mbox{for any}\ t\in T_{i}^{+}\right\\},$
$\displaystyle\Lambda_{i1}^{-}=\left\\{\gamma\bigg{|}\frac{\mathrm{d}C_{it}(\gamma)}{\mathrm{d}t}\geq\Phi^{-1}(1-\alpha)\
\mbox{for any}\ t\in T_{i}^{-}\right\\},$

$i=1,2,\cdots,m$, where

$\Phi^{-1}(\alpha)=\frac{\sqrt{3}}{\pi}\mathrm{ln}\frac{\alpha}{1-\alpha}.$

Because $T_{i}^{+}$ and $T_{i}^{-}$ are disjoint sets, and
$C_{1t},\cdots,C_{mt}$ are independent increment processes, we have

$\mbox{\euler M}\\{\Lambda_{i1}^{+}\\}=\alpha,\quad\mbox{\euler
M}\\{\Lambda_{i1}^{-}\\}=\alpha,\quad\mbox{\euler
M}\\{\Lambda_{i1}^{+}\cap\Lambda_{i1}^{-}\\}=\alpha,$

$i=1,2,\cdots,m$. For any $\gamma\in\Lambda_{i1}^{+}\cap\Lambda_{i1}^{-}$, it
is apparent that for any $t$,

$\displaystyle
g_{i}\left(t,X_{t}(\gamma),\cdots,\frac{\mathrm{d}^{n-1}X_{t}(\gamma)}{\mathrm{d}t^{n-1}}\right)\frac{\mathrm{d}C_{it}(\gamma)}{\mathrm{d}t}$
$\displaystyle\leq\left|g_{i}\left(t,X_{t}^{\alpha},\cdots,\frac{\mathrm{d}^{n-1}X_{t}^{\alpha}}{\mathrm{d}t^{n-1}}\right)\right|\Phi^{-1}(\alpha),$

$i=1,2,\cdots,m$. Let
$\Lambda_{1}^{+}\cap\Lambda_{1}^{-}=\bigcap\limits_{i=1}^{m}(\Lambda_{i1}^{+}\cap\Lambda_{i1}^{-})$,
$i=1,2,\cdots,m$.

Since $C_{1t},C_{2t},\cdots,C_{mt}$ are independent and $\mbox{\euler
M}\\{\Lambda_{i1}^{+}\cap\Lambda_{i1}^{-}\\}=\alpha$, $i=1,2,\cdots,m$, we
have

$\mbox{\euler M}\\{\Lambda_{1}^{+}\cap\Lambda_{1}^{-}\\}=\mbox{\euler
M}\left\\{\bigcap\limits_{i=1}^{m}(\Lambda_{i1}^{+}\cap\Lambda_{i1}^{-})\right\\}=\bigwedge_{1\leq
i\leq m}\mbox{\euler
M}\big{\\{}\Lambda_{i1}^{+}\cap\Lambda_{i1}^{-}\big{\\}}=\alpha.$

So, for any $\gamma\in\Lambda_{1}^{+}\cap\Lambda_{1}^{-}$, we get for any $t$,

$\displaystyle\sum_{i=1}^{m}g_{i}\left(\\!t,X_{t}(\\!\gamma),\cdots,\frac{\mathrm{d}^{n\\!-\\!1}X_{t}(\\!\gamma)}{\mathrm{d}t^{n\\!-\\!1}}\\!\right)\frac{\mathrm{d}C_{it}(\gamma)}{\mathrm{d}t}$
$\displaystyle\leq\sum_{i=1}^{m}\left|g_{i}\left(\\!t,X_{t}^{\alpha},\cdots,\frac{\mathrm{d}^{n\\!-\\!1}X_{t}^{\alpha}}{\mathrm{d}t^{n\\!-\\!1}}\\!\right)\right|\Phi^{\\!-\\!1}(\alpha),$

i.e.

$\displaystyle
f\left(t,X_{t}(\gamma),\cdots,\frac{\mathrm{d}^{n-1}X_{t}(\gamma)}{\mathrm{d}t^{n-1}}\right)+\sum_{i=1}^{m}g_{i}\left(t,X_{t}(\gamma),\cdots,\frac{\mathrm{d}^{n-1}X_{t}(\gamma)}{\mathrm{d}t^{n-1}}\right)\frac{\mathrm{d}C_{it}(\gamma)}{\mathrm{d}t}$
$\displaystyle\leq
f\left(t,X_{t}^{\alpha},\cdots,\frac{\mathrm{d}^{n-1}X_{t}^{\alpha}}{\mathrm{d}t^{n-1}}\right)+\sum_{i=1}^{m}\left|g_{i}\left(t,X_{t}^{\alpha},\cdots,\frac{\mathrm{d}^{n-1}X_{t}^{\alpha}}{\mathrm{d}t^{n-1}}\right)\right|\Phi^{-1}(\alpha).$

Since

$f\left(t,X_{t}^{\alpha},\cdots,\frac{\mathrm{d}^{n-1}X_{t}^{\alpha}}{\mathrm{d}t^{n-1}}\right)+\sum\limits_{i=1}\limits^{m}\left|g_{i}\left(t,X_{t}^{\alpha},\cdots,\frac{\mathrm{d}^{n-1}X_{t}^{\alpha}}{\mathrm{d}t^{n-1}}\right)\right|\Phi^{-1}(\alpha)$

is a monotonically increasing function with respect to

$X_{t}^{\alpha},\frac{\mathrm{d}X_{t}^{\alpha}}{\mathrm{d}t},\cdots,\frac{\mathrm{d}^{n-2}X_{t}^{\alpha}}{\mathrm{d}t^{n-2}},$

according to Theorem 2.2, we have

$X_{t}\leq X_{t}^{\alpha},\ \forall t.$

Note that $\Lambda_{1}^{+}\cap\Lambda_{1}^{-}\subset\\{X_{t}\leq
X_{t}^{\alpha},\forall t\\}$, we can get

$\mbox{\euler M}\\{X_{t}\leq X_{t}^{\alpha},\forall t\\}\geq\mbox{\euler
M}\\{\Lambda_{1}^{+}\cap\Lambda_{1}^{-}\\}=\alpha.$ (17)

Next, let

$\displaystyle\Lambda_{i2}^{+}=\left\\{\gamma\bigg{|}\frac{\mathrm{d}C_{it}(\gamma)}{\mathrm{d}t}>\Phi^{-1}(\alpha)\
\mbox{for any}\ t\in T_{i}^{+}\right\\},$
$\displaystyle\Lambda_{i2}^{-}=\left\\{\gamma\bigg{|}\frac{\mathrm{d}C_{it}(\gamma)}{\mathrm{d}t}<\Phi^{-1}(1-\alpha)\
\mbox{for any}\ t\in T_{i}^{-}\right\\},$

$i=1,2,\cdots,m$. Because $T_{i}^{+}$ and $T_{i}^{-}$ are disjoint sets and
$C_{1t},\cdots,C_{mt}$ are independent increment processes, we get

$\mbox{\euler M}\\{\Lambda_{i2}^{+}\\}=1-\alpha,\quad\mbox{\euler
M}\\{\Lambda_{i2}^{-}\\}=1-\alpha,\quad\mbox{\euler
M}\\{\Lambda_{i2}^{+}\cap\Lambda_{i2}^{-}\\}=1-\alpha,$

$i=1,2,\cdots,m$. Considering
$\forall\gamma\in\Lambda_{i2}^{+}\cap\Lambda_{i2}^{-}$, it is apparent that
for any $t$,

$\displaystyle
g_{i}\left(t,X_{t}(\gamma),\cdots,\frac{\mathrm{d}^{n-1}X_{t}(\gamma)}{\mathrm{d}t^{n-1}}\right)\frac{\mathrm{d}C_{it}(\gamma)}{\mathrm{d}t}$
$\displaystyle>\left|g_{i}\left(t,X_{t}^{\alpha},\cdots,\frac{\mathrm{d}^{n-1}X_{t}^{\alpha}}{\mathrm{d}t^{n-1}}\right)\right|\Phi^{-1}(\alpha),$

$i=1,2,\cdots,m$.

Let
$\Lambda_{2}^{+}\cap\Lambda_{2}^{-}=\bigcap\limits_{i=1}^{m}(\Lambda_{i2}^{+}\cap\Lambda_{i2}^{-}),i=1,2,\cdots,m$.
Since $C_{1t},C_{2t},\cdots,C_{mt}$ are independent and $\mbox{\euler
M}\\{\Lambda_{i2}^{+}\cap\Lambda_{i2}^{-}\\}\\!=\\!1-\alpha,i=1,2,\\!\cdots\\!,m$,
we have

$\mbox{\euler M}\\{\Lambda_{2}^{+}\cap\Lambda_{2}^{-}\\}=\mbox{\euler
M}\left\\{\bigcap\limits_{i=1}^{m}(\Lambda_{i2}^{+}\cap\Lambda_{i2}^{-})\right\\}=\bigwedge_{1\leq
i\leq m}\mbox{\euler
M}\big{\\{}\Lambda_{i2}^{+}\cap\Lambda_{i2}^{-}\big{\\}}=1-\alpha.$

So, for any $\gamma\in\Lambda_{2}^{+}\cap\Lambda_{2}^{-}$, we get for any $t$,

$\displaystyle\sum_{i=1}^{m}g_{i}\left(\\!t,X_{t}(\\!\gamma),\cdots,\frac{\mathrm{d}^{n\\!-\\!1}X_{t}(\\!\gamma)}{\mathrm{d}t^{n\\!-\\!1}}\\!\right)\frac{\mathrm{d}C_{it}(\gamma)}{\mathrm{d}t}$
$\displaystyle>\sum_{i=1}^{m}\left|g_{i}\left(\\!t,X_{t}^{\alpha},\cdots,\frac{\mathrm{d}^{n\\!-\\!1}X_{t}^{\alpha}}{\mathrm{d}t^{n\\!-\\!1}}\\!\right)\right|\Phi^{\\!-\\!1}(\alpha).$

According to Theorem 2.1, we have

$X_{t}>X_{t}^{\alpha},\ \forall t.$

It is obvious that
$\Lambda_{2}^{+}\cap\Lambda_{2}^{-}\subset\\{X_{t}>X_{t}^{\alpha},\forall
t\\}$. Hence

$\mbox{\euler M}\\{X_{t}>X_{t}^{\alpha},\forall t\\}\geq\mbox{\euler
M}\\{\Lambda_{2}^{+}\cap\Lambda_{2}^{-}\\}=1-\alpha.$ (18)

Since the opposite of $\\{X_{t}\leq X_{t}^{\alpha},\forall t\\}$ is
$\\{X_{t}\nleq X_{t}^{\alpha},\forall t\\}$, on the basis of duality axiom, we
can get

$\mbox{\euler M}\\{X_{t}\leq X_{t}^{\alpha},\forall t\\}+\mbox{\euler
M}\\{X_{t}\nleq X_{t}^{\alpha},\forall t\\}=1.$

Besides, $\\{X_{t}>X_{t}^{\alpha},\forall t\\}\subset\\{X_{t}\nleq
X_{t}^{\alpha},\forall t\\}$ implies

$\mbox{\euler M}\\{X_{t}\leq X_{t}^{\alpha},\forall t\\}+\mbox{\euler
M}\\{X_{t}>X_{t}^{\alpha},\forall t\\}\leq 1.$ (19)

Thus, from the (17),(18) and (19), it is evident that

$\displaystyle\mbox{\euler M}\\{X_{t}\leq X_{t}^{\alpha},\forall t\\}=\alpha,$
$\displaystyle\mbox{\euler M}\\{X_{t}>X_{t}^{\alpha},\forall t\\}=1-\alpha.$

The theorem is proved.

Through the discussion of the above theorems, we can easily draw the following
theorem.

######  Theorem 2.4.

Let $f\big{(}t,z_{0},z_{1},\cdots,z_{n-1}\big{)}$ be a function on
$D=[0,a]\\!\times\\!\mathbb{D}_{0}\subset\mathbb{R}^{n\\!+\\!1}$
$(\mbox{where}\ \mathbb{D}_{0}\ \mbox{is an n-dimensional bounded region})$
satisfying local Lipschitz conditions in $z_{0},z_{1},\cdots,z_{n-1}$, and
$g_{1}(t),g_{2}(t),\dots,g_{m}(t)$ be integrable functions with respect to
$t$. Assume that $f\big{(}t,z_{0},z_{1},\cdots,z_{n-1}\big{)}$ is a
monotonically increasing function with respect to $z_{0},z_{1},\cdots,$
$z_{n-2}$. If $X_{t}$ and $X_{t}^{\alpha}$ are the solution of

$\frac{\mathrm{d}^{n}X_{t}}{\mathrm{d}t^{n}}=f\bigg{(}t,X_{t},\cdots,\frac{\mathrm{d}^{n-1}X_{t}}{\mathrm{d}t^{n-1}}\bigg{)}+\sum_{i=1}^{m}g_{i}(t)\frac{\mathrm{d}C_{it}}{\mathrm{d}t}$
(20)

and

$\frac{\mathrm{d}^{n}X_{t}^{\alpha}}{\mathrm{d}t^{n}}=f\bigg{(}t,X_{t}^{\alpha},\cdots,\frac{\mathrm{d}^{n-1}X_{t}^{\alpha}}{\mathrm{d}t^{n-1}}\bigg{)}+\sum_{i=1}^{m}|g_{i}(t)|\Phi^{-1}(\alpha),$
(21)

respectively, where

$\Phi^{-1}(\alpha)=\frac{\sqrt{3}}{\pi}\mathrm{ln}\frac{\alpha}{1-\alpha},\qquad
0<\alpha<1,$

then $X_{t}$ has an inverse uncertainty distribution

$\Psi_{t}^{-1}(\alpha)=X_{t}^{\alpha}.$ (22)

######  Example 2.2.

Suppose that $X_{t}$ and $X_{t}^{\alpha}$ are the solution of

$\left\\{\begin{aligned}
&\frac{\mathrm{d}^{2}X_{t}}{\mathrm{d}t^{2}}=2\frac{\mathrm{d}X_{t}}{\mathrm{d}t}+3X_{t}+e^{-t}\frac{\mathrm{d}C_{t}}{\mathrm{d}t}\\\
&X_{(0)}=0\\\ &X^{\prime}_{(0)}=0\end{aligned}\right.$ (23)

and

$\left\\{\begin{aligned}
&\frac{\mathrm{d}^{2}X^{\alpha}_{t}}{\mathrm{d}t^{2}}=2\frac{\mathrm{d}X^{\alpha}_{t}}{\mathrm{d}t}+3X^{\alpha}_{t}+e^{-t}\Phi^{-1}(\alpha)\\\
&X^{\alpha}_{(0)}=0\\\ &{X^{\alpha}_{(0)}}^{\prime}=0\end{aligned}\right.$
(24)

where

$\Phi^{-1}(\alpha)=\frac{\sqrt{3}}{\pi}\mathrm{ln}\frac{\alpha}{1-\alpha},\qquad
0<\alpha<1,$

respectively.

Note that

$2\frac{\mathrm{d}X^{\alpha}_{t}}{\mathrm{d}t}+3X^{\alpha}_{t}+e^{-t}\Phi^{-1}(\alpha)$

is a monotonically increasing function with respect to $X_{t}^{\alpha}$, so
$X^{\alpha}_{t}$ is the $\alpha$-path of $X_{t}$. Solve this Eq.(24) through
Euler exponential function method and obtain the solution

$X_{t}^{\alpha}=\frac{\sqrt{3}}{16\pi}\bigg{(}e^{3t}-e^{-t}-4te^{-t}\bigg{)}\cdot\ln\frac{\alpha}{1-\alpha},\quad
0<\alpha<1.$

Figure 3: Different values of $\alpha$ of $\Psi_{t}^{-1}(\alpha)$

According to Theorem 2.4, $X_{t}^{\alpha}$ is the $\alpha$-path of $X_{t}$
(shown in Fig.3), and the inverse uncertainty distribution of $X_{t}$ is

$\Psi_{t}^{-1}(\alpha)=X_{t}^{\alpha}.$

##  3 Residual

In general, the observations of an uncertain process are discrete points when
the time intervals are not very short. In order to make a connection between
HUDE and observations, we first define the residual of HUDE. Let us consider a
high-order uncertain differential equation

$\frac{\mathrm{d}^{n}X_{t}}{\mathrm{d}t^{n}}=f\bigg{(}t,X_{t},\frac{\mathrm{d}X_{t}}{\mathrm{d}t},\cdots,\frac{\mathrm{d}^{n-1}X_{t}}{\mathrm{d}t^{n-1}}\bigg{)}+\sum_{i=1}^{m}g_{i}\bigg{(}t,X_{t},\frac{\mathrm{d}X_{t}}{\mathrm{d}t},\cdots,\frac{\mathrm{d}^{n-1}X_{t}}{\mathrm{d}t^{n-1}}\bigg{)}\frac{\mathrm{d}C_{it}}{\mathrm{d}t}$

where $f$ and $g_{i}\ (i=1,2,\cdots,m)$ are known continuous functions and
$C_{1t},C_{2t},\cdots,$ $C_{mt}$ are independent Liu processes. Assume that

$x_{t_{j}},\ x^{\prime}_{t_{j}},\ ,\cdots,\ x^{(n\\!-\\!1)}_{t_{j}}$ (25)

are observations of
$X_{t},\frac{\mathrm{d}X_{t}}{\mathrm{d}t},\cdots,\frac{\mathrm{d}^{n-1}X_{t}}{\mathrm{d}t^{n-1}}$
at $t_{j}\ (j=1,2,\cdots,l)$, where $t_{1}<t_{2}<\cdots<t_{l}$. Observation
$x_{t_{l+1}}$ at time $t_{l+1}$ is obtained. For convenience, all observed
data is listed in Table 1.

For any given index $j$ with $1\leq j\leq l$, we consider solving the updated
high-order uncertain differential equation,

$\left\\{\begin{array}[]{l}\vspace{1ex}\displaystyle\frac{\mathrm{d}^{n}X_{t}}{\mathrm{d}t^{n}}\\!\\!=\\!\\!f\bigg{(}t,X_{t},\frac{\mathrm{d}X_{t}}{\mathrm{d}t},\cdots,\frac{\mathrm{d}^{n\\!-\\!1}X_{t}}{\mathrm{d}t^{n\\!-\\!1}}\bigg{)}\\!+\\!\sum\limits_{i=1}^{m}g_{i}\bigg{(}t,X_{t},\frac{\mathrm{d}X_{t}}{\mathrm{d}t},\cdots,\frac{\mathrm{d}^{n\\!-\\!1}X_{t}}{\mathrm{d}t^{n\\!-\\!1}}\bigg{)}\frac{\mathrm{d}C_{it}}{\mathrm{d}t}\\\
\vspace{1ex}X_{t_{j}}=x_{t_{j}}\\\
\displaystyle\frac{\mathrm{d}X_{t}}{\mathrm{d}t}\bigg{|}_{t=t_{j}}=x^{\prime}_{t_{j}}\\\
\qquad\quad\vdots\\\
\displaystyle\frac{\mathrm{d}^{n-1}X_{t}}{\mathrm{d}t^{n-1}}\bigg{|}_{t=t_{j}}=x^{(n\\!-\\!1)}_{t_{j}}\\\
\end{array}\right.$ (26)

where $x_{t_{j}},x^{\prime}_{t_{j}},\cdots,x^{(n\\!-\\!1)}_{t_{j}}$ are
observations at time $t_{j}$. The uncertainty distribution of $X_{t_{j+1}}$
initialized at $t_{j}$, denoted by $\Phi_{t_{j+1}}$, could be obtained by
solving Eq.(26).

Table 1: Observed data at time $t_{j}$ ($1\leq j\leq l+1$) $t_{1}$ | $t_{2}$ | $\cdots$ | $t_{l}$ | $t_{l+1}$
---|---|---|---|---
$x_{t_{1}}$ | $x_{t_{2}}$ | $\cdots$ | $x_{t_{l}}$ | $x_{t_{l+1}}$
$x^{\prime}_{t_{1}}$ | $x^{\prime}_{t_{2}}$ | $\cdots$ | $x^{\prime}_{t_{l}}$ |
$\vdots$ | $\vdots$ | $\vdots$ | $\vdots$ |
$x_{t_{1}}^{(n\\!-\\!1)}$ | $x^{(n\\!-\\!1)}_{t_{2}}$ | $\cdots$ | $x^{(n\\!-\\!1)}_{t_{l}}$ |

######  Proposition 3.1.

Let $\xi$ be a uncertain variable with regular uncertainty distribution
$\Phi(x)$, then $\Phi(\xi)$ is a linear uncertain variable $\mathcal{L}(0,1)$.

For $1\leq j\leq l$, we know the uncertain distribution of $X_{t_{j+1}}$ and
its observed value $x_{t_{j+1}}$. Then the $j$th residual is defined as

$\varepsilon_{j}=\Phi_{t_{j+1}}(x_{t_{j+1}}).$ (27)

Then we have a total of $l$ residuals

$\varepsilon_{1},~{}\varepsilon_{2},~{}\cdots,~{}\varepsilon_{l}.$

By Proposition 3.1, all the residuals form a sample of linear uncertainty
distribution $\mathcal{L}(0,1)$.

For most of high-order uncertain differential equations, it is not easy to
obtain analytical solutions. Sometimes, it is even unpractical to calculate
any analytical solutions. In this case, we need to use numerical methods to
calculate the residual. The following algorithm is able to calculate the $j$th
residual $\varepsilon_{j}$.

Algorithm 1: Numerical method for calculating residuals
Step 0: Set $l=0$, $r=1$ and a precision $\delta=0.0001$.
Step 1: Set $\displaystyle\alpha=(l+r)/2$.
Step 2: Compute $X_{t_{j+1}}^{\alpha}$ of the updated Eq.(26) by Euler method.
Step 3: If $X_{t_{j+1}}^{\alpha}<x_{t_{j+1}}$, then $l=\alpha$. Otherwise,
$r=\alpha$.
Step 4: If $|l-r|>\delta$, then go to Step 1.
Step 5: Output $\varepsilon_{j}=(l+r)/2$.

######  Remark 3.1.

However, in general, we can only obtain the observations of $X_{t}$ at
different times, and it is difficult to get the value of
$X^{\prime}_{t},\cdots,X^{(n-1)}_{t}$ at the corresponding time by
observation. In this case, suppose

$x_{t_{1}},x_{t_{2}},\cdots,x_{t_{l+n-1}}$ (28)

are $l+n-1$ observations of $X_{t}$ at the times
$t_{1},t_{2},\cdots,t_{l+n-1}$ with $t_{1}<t_{2}<\cdots<t_{l+n-1}$,
respectively. Define

$x^{(v)}_{t_{j}}=\frac{x^{(v-1)}_{t_{j+1}}-x^{(v-1)}_{t_{j}}}{t_{j+1}-t_{j}}$
(29)

is the value of $X_{t}^{(v)}(v=1,2,\cdots,n-1)$ at the $t_{j}$ where $1\leq
j\leq l+n-2$. Therefore, $l+n-1$ data are able to get the complete initial
information of the first $l$ moments. All data is shown in Table 2.

Table 2: Observed data at time $t_{j}$ ($1\leq j\leq l+n-1$) $t_{1}$ | $t_{2}$ | $\cdots$ | $t_{l}$ | $t_{l+1}$ | $\cdots$ | $t_{l+n}$ | $t_{l+n-1}$
---|---|---|---|---|---|---|---
$x_{t_{1}}$ | $x_{t_{2}}$ | $\cdots$ | $x_{t_{l}}$ | $x_{t_{l+1}}$ | $\cdots$ | $x_{t_{l+n}}$ | $x_{t_{l+n-1}}$
$x^{\prime}_{t_{1}}$ | $x^{\prime}_{t_{2}}$ | $\cdots$ | $x^{\prime}_{t_{l}}$ | $x^{\prime}_{t_{l+1}}$ | $\cdots$ | $x^{\prime}_{t_{l+n}}$ |
$\vdots$ | $\vdots$ |  | $\vdots$ | $\vdots$ |  |  |
$x_{t_{1}}^{(n\\!-\\!2)}$ | $x^{(n\\!-\\!2)}_{t_{2}}$ | $\cdots$ | $x^{(n\\!-\\!2)}_{t_{l}}$ | $x^{(n\\!-\\!2)}_{t_{l+1}}$ |  |  |
$x_{t_{1}}^{(n\\!-\\!1)}$ | $x^{(n\\!-\\!1)}_{t_{2}}$ | $\cdots$ | $x^{(n\\!-\\!1)}_{t_{l}}$ |  |  |  |

Using difference (29) is the simplest form, in addition, derivative can be
approximated by numerical differentiation or Lagrange interpolation process.
For example, consider central difference

$x^{(v)}_{t_{j}}=\frac{x^{(v-1)}_{t_{j+1}}-x^{(v-1)}_{t_{j-1}}}{t_{j+1}-t_{j-1}}.$

There are two main use of residuals: parameter estimation and hypothesis test,
which will be discussed in the next two sections.

##  4 Parameter estimation

When we apply the model of HUDE to practical problems, the differential
equations normally contain unknown parameters. In this section, we will
estimate unknown parameters based on observed data and residuals.

Consider the following high-order uncertain differential equation

$\frac{\mathrm{d}^{n}X_{t}}{\mathrm{d}t^{n}}=f\bigg{(}t,X_{t},\cdots,\frac{\mathrm{d}^{n\\!-\\!1}X_{t}}{\mathrm{d}t^{n\\!-\\!1}};\bm{\theta}\bigg{)}+\sum_{i=1}^{m}g_{i}\bigg{(}t,X_{t},\cdots,\frac{\mathrm{d}^{n\\!-\\!1}X_{t}}{\mathrm{d}t^{n\\!-\\!1}};\bm{\theta}\bigg{)}\frac{\mathrm{d}C_{it}}{\mathrm{d}t}$
(30)

where $f$ and $g_{i}\ (i=1,2,\cdots,m)$ are known continuous functions,
$C_{1t},C_{2t},\cdots,C_{mt}$ are independent Liu processes, and
$\bm{\theta}(\bm{\theta}\in\mathcal{I},\mathcal{I}\subseteq\mathbb{R}^{p})$ is
an unknown $p$-vector of parameters. Suppose the observed data is listed in
Table 1.

Then for every $j$ with $1\leq j\leq l$, the updated high-order uncertain
differential equation

$\left\\{\begin{array}[]{l}\vspace{1ex}\displaystyle\frac{\mathrm{d}^{n}X_{t}}{\mathrm{d}t^{n}}\\!\\!=\\!\\!f\bigg{(}t,X_{t},\frac{\mathrm{d}X_{t}}{\mathrm{d}t},\\!\cdots\\!,\frac{\mathrm{d}^{n\\!-\\!1}X_{t}}{\mathrm{d}t^{n\\!-\\!1}};\\!\bm{\theta}\bigg{)}\\!+\\!\sum\limits_{i=1}^{m}g_{i}\bigg{(}t,X_{t},\frac{\mathrm{d}X_{t}}{\mathrm{d}t},\\!\cdots\\!,\frac{\mathrm{d}^{n\\!-\\!1}X_{t}}{\mathrm{d}t^{n\\!-\\!1}};\\!\bm{\theta}\bigg{)}\frac{\mathrm{d}C_{it}}{\mathrm{d}t}\\\
\vspace{1ex}X_{t_{j}}=x_{t_{j}}\vspace{1ex}\\\
\displaystyle\frac{\mathrm{d}X_{t}}{\mathrm{d}t}\bigg{|}_{t=t_{j}}=x^{\prime}_{t_{j}}\\\
\qquad\quad\vdots\\\
\displaystyle\frac{\mathrm{d}^{n-1}X_{t}}{\mathrm{d}t^{n-1}}\bigg{|}_{t=t_{j}}=x^{(n\\!-\\!1)}_{t_{j}}\\\
\end{array}\right.$ (31)

contains unknown vector $\bm{\theta}$. Then by Eq. (27), the $j$th residual is

$\varepsilon_{j}(\bm{\theta})=\Phi_{t_{j+1}}(x_{t_{j+1}}|\bm{\theta}).$ (32)

By Proposition 3.1,

$\varepsilon_{1}(\bm{\theta}),\varepsilon_{2}(\bm{\theta}),\\!\cdots\\!,\varepsilon_{l}(\bm{\theta})$

form a sample of $\mathcal{L}(0,1)$. Based on the residuals, there are two
methods to estimate the unknown $p$-vector $\bm{\theta}$.

The first esitmation of $\bm{\theta}$ is the moment estimation. As the $k$th
population moment of the linear uncertainty distribution $\mathcal{L}(0,1)$ is

$\frac{1}{k+1}.$

According to the principle that the $k$th sample moment is equal to the $k$th
population moment, the moment estimate $\bm{\theta}$ should resolve the system
of equations

$\frac{1}{l}\sum_{j=1}^{l}\varepsilon^{k}_{j}(\bm{\theta})=\frac{1}{k+1},\quad
k=1,2,\cdots,p,$ (33)

where $p$ is the dimension of $\bm{\theta}$. Thus, the solution of the
following minimization problem,

$\left\\{\begin{array}[]{ll}\vspace{1mm}\min\limits_{\bm{\theta}}\sum\limits_{k=1}^{p}\bigg{(}\displaystyle\frac{1}{l}\sum\limits_{j=1}^{l}\varepsilon^{k}_{j}(\bm{\theta})-\frac{1}{k+1}\bigg{)}^{2}\\\
\rm{subject\ to:}\\\ \
\qquad\qquad\bm{\theta}\in\mathcal{I}\end{array}\right.$ (34)

is the moment estimation of $\bm{\theta}$. This minimization problem (34)
could be solved by MATLAB111MATLAB R2021a, 9.10.0.1602886, maci64,
Optimization Toolbox, “fminsearch” function..

######  Remark 4.1.

The optimal value of the objective function in Eq (34) should be very close to
zero. In the actual calculation using MATLAB, it is generally considered that
the value of objective function should be less than $10^{-10}$. Otherwise, we
can assume that data do not fit the proposed uncertain differential equation.

The other estimation is the maximum likelihood estimation proposed by [9].
Given a detection level $\alpha$, the maximum likelihood estimation of
$\bm{\theta}$ is the solution of the following system of equations,

$\begin{cases}\displaystyle\varepsilon^{\prime}_{i^{*}(\bm{\theta})}(\bm{\theta})=\frac{\alpha}{2}&\\\
\displaystyle\varepsilon^{\prime}_{i^{*}(\bm{\theta})+\lceil
l(1-\alpha)\rceil-1}(\bm{\theta})=1-\frac{\alpha}{2}&\\\
i^{*}(\bm{\theta})=\mathop{\arg\min}\limits_{1\leq i\leq l-\lceil
l(1-\alpha)\rceil+2}\varepsilon^{\prime}_{i+\lceil
l(1-\alpha)\rceil-1}(\bm{\theta})-\varepsilon^{\prime}_{i}(\bm{\theta}),\end{cases}$

where

$\\{\varepsilon^{\prime}_{1}(\bm{\theta}),\varepsilon^{\prime}_{2}(\bm{\theta}),\cdots,\varepsilon^{\prime}_{l}(\bm{\theta})\\}$

with

$\varepsilon^{\prime}_{1}(\bm{\theta})\leq\varepsilon^{\prime}_{2}(\bm{\theta}),\cdots\leq\varepsilon^{\prime}_{l}(\bm{\theta})$

is a rearrangement of

$\\{\varepsilon_{1}(\bm{\theta}),\varepsilon_{2}(\bm{\theta}),\cdots,\varepsilon_{l}(\bm{\theta})\\}.$

##  5 Uncertain hypothesis test

After the unknown parameters have been estimated using the method in Sect.4,
it is crucial to test whether a HUDE is a good fit to the observed datas. Here
we will use the uncertain hypothesis test introdueced by [19] to evaluate the
fitness of HUDE.

Consider an high-order uncertain differential equation

$\frac{\mathrm{d}^{n}X_{t}}{\mathrm{d}t^{n}}=f\bigg{(}t,X_{t},\cdots,\frac{\mathrm{d}^{n\\!-\\!1}X_{t}}{\mathrm{d}t^{n\\!-\\!1}};\bm{\theta}\bigg{)}+\sum_{i=1}^{m}g_{i}\bigg{(}t,X_{t},\cdots,\frac{\mathrm{d}^{n\\!-\\!1}X_{t}}{\mathrm{d}t^{n\\!-\\!1}};\bm{\theta}\bigg{)}\frac{\mathrm{d}C_{it}}{\mathrm{d}t}$
(35)

where $f$ and $g_{i}\ (i=1,2,\cdots,m)$ are known continuous functions and
$C_{1t},C_{2t},\cdots,$ $C_{mt}$ are independent Liu processes but
$\bm{\theta}(\bm{\theta}\in\mathcal{I},\mathcal{I}\subseteq\mathbb{R}^{p})$ is
an unknown $p$-vector of parameters. Assume that

$x_{t_{1}},x_{t_{2}},\cdots,x_{t_{l}}$ (36)

are $l$ observations of $X_{t}$ at the times $t_{1},t_{2},\cdots,t_{l}$ with
$t_{1}<t_{2}<\cdots<t_{l}$, where $l>n$, respectively. For each index $j(1\leq
j\leq l-n+1)$, we solve the updated high-order uncertain differential equation

$\left\\{\begin{array}[]{l}\vspace{1ex}\displaystyle\frac{\mathrm{d}^{n}X_{t}}{\mathrm{d}t^{n}}\\!=\\!f\bigg{(}t,X_{t},\frac{\mathrm{d}X_{t}}{\mathrm{d}t},\\!\cdots\\!,\frac{\mathrm{d}^{n-1}X_{t}}{\mathrm{d}t^{n-1}};\bm{\theta}\bigg{)}\\!+\\!\sum\limits_{i=1}^{m}g_{i}\bigg{(}t,X_{t},\frac{\mathrm{d}X_{t}}{\mathrm{d}t},\\!\cdots\\!,\frac{\mathrm{d}^{n\\!-\\!1}X_{t}}{\mathrm{d}t^{n\\!-\\!1}};\bm{\theta}\bigg{)}\frac{\mathrm{d}C_{it}}{\mathrm{d}t}\\\
\vspace{1ex}X_{t_{j}}=x_{t_{j}}\\\
\displaystyle\frac{\mathrm{d}X_{t}}{\mathrm{d}t}\bigg{|}_{t=t_{j}}=x^{\prime}_{t_{j}}\\\
\qquad\quad\vdots\\\
\displaystyle\frac{\mathrm{d}^{n-1}X_{t}}{\mathrm{d}t^{n-1}}\bigg{|}_{t=t_{j}}=x^{(n-1)}_{t_{j}}\\\
\end{array}\right.$ (37)

where $x_{t_{j}},x^{\prime}_{t_{j}},\cdots,x^{(n\\!-\\!1)}_{t_{j}}$ are $n$
new initial values according to the form of difference (29) at the new initial
time $t_{j}$ with $1\leq j\leq l-n+1$, respectively.

For any given $\bm{\theta}$, on the basis of Sec. 4, we can produce $l-n+1$
residuals

$\varepsilon_{1}(\bm{\theta}),\varepsilon_{2}(\bm{\theta}),\\!\cdots\\!,\varepsilon_{l\\!-\\!n+1}(\bm{\theta})$

of Eq.(37) corresponding to the observed data (36).

If the high-order uncertain differential Eq.(35) does fit the observed data
(36), then

$\varepsilon_{1},\varepsilon_{2},\\!\cdots\\!,\varepsilon_{l\\!-\\!n+1}\sim\mathcal{L}(0,1).$

That is, testing whether a high-order uncertain differential equation fits the
observed data well is equivalent to testing whether these residuals
$\varepsilon_{1},\varepsilon_{2},\\!\cdots\\!,\varepsilon_{l\\!-\\!n+1}$ fit
the uncertainty distribution $\mathcal{L}(0,1)$, i.e.

$\varepsilon_{1},\varepsilon_{2},\\!\cdots\\!,\varepsilon_{l\\!-\\!n+1}\sim\mathcal{L}(0,1).$

Given a significance level $\alpha\ (\text{e.g.}\ 0.05)$, the test is

$\displaystyle
W=\bigg{\\{}(\varepsilon_{1},\varepsilon_{2},\cdots,\varepsilon_{l\\!-\\!n+1}):$
$\displaystyle{\text{there are at least $\alpha$ of indexes j's with }}1\leq
j\leq l\\!-\\!n\\!+\\!1$ $\displaystyle{\text{ such that
}}\varepsilon_{j}<\frac{\alpha}{2}{\text{ or
}}\varepsilon_{j}>1-\frac{\alpha}{2}\bigg{\\}}.$

If the vector of the $l-n+1$ residuals
$\varepsilon_{1},\varepsilon_{2},\cdots,\varepsilon_{l\\!-\\!n+1}$ belongs to
the test $W$, i.e.,

$(\varepsilon_{1},\varepsilon_{2},\cdots,\varepsilon_{l\\!-\\!n+1})\in W,$

then Eq.(35) is not a good fit to the observed data (36).

If

$(\varepsilon_{1},\varepsilon_{2},\cdots,\varepsilon_{l\\!-\\!n+1})\notin W,$

then Eq.(35) is a good fit to the observed data (36).

##  6 An application in nuclear reactor kinetics

Nuclear safety has emerged as a crucial consideration for many nations in
their pursuit of peaceful development. Despite the universal commitment to
nuclear disarmament and non-proliferation, research in the nuclear industry
continues unabated.

The nuclear reactor kinetics equations have been studied and modeled by many
scholars([2]; [1, 3]). The nuclear reactor kinetics is also known as neutron
population kinetics which studies the dynamic change of neutron population in
reactor. Neutron population determines the change of power level over time and
are affected by the control rod position and other factors([2]).

Due to the specific operational complexities of nuclear reactors, research
into refining nuclear reactor models is imperative. Furthermore, in cases
involving newly discovered radioactive isotopes or neutrons lacking
comprehensive data, experts often rely on past experience to give credibility
to analyze their properties. The current nuclear reactor models do not
adequately account for such uncertainties.

###  6.1 Nuclear reactor kinetics driven by Liu process

In this part, the quantitative relationship between neutron population and
time variation in a nuclear reactor is modeled. In a real nuclear reactor, if
the delayed neutrons are not taken into account, the neutron population
increases so rapidly in the supercritical state that the reactor cannot be
controlled, which is extremely dangerous for a nuclear reactor. A delayed
neutron was produced by a delayed neutron precursor. Every delayed neutron
precursor likes Br-87, Uranium-235 is produced at the instant of fission and
releases delayed neutrons after a slight delay. In order to better establish
the nuclear reactor kinetics model, the nuclear reactor is supposed to be
large enough that energy and space effects can be ignored.

Figure 4: Dynamical system of a nuclear reactor

Normally, there is more than one delayed neutron precursors in a reactor, and
nuclear physicists divide these delayed neutron precursors into roughly six
groups based on average lifetime and decay constant. Therefore, it is assumed
that there are six groups of delayed neutron precursors in this reactor. Fig.4
depicts the dynamic processes of the neutron concentrations and precursors
denisty.

In this figure, an extraneous neutron source $B$ transports the neutron at a
constant rate to the nuclear reactor. $N(t)$ is the neutron population with
respect to time $t$ and the variable $Q_{i}(t)$ is $i$-th delayed neutron
precursor density over time. The following parameters are important to
describe the dynamical system. Assume that $k$ is effective multiplication
constant of neutron which is calculated as neutrons population in the new
generation divided by neutrons population in the old generation and neutron
lifetime is $l$. Consider that $\beta_{i}$ is the fraction of $i$-th delayed
neutron precursor and $\beta=\sum\limits_{i=1}\limits^{6}\beta_{i}$ is the
delayed neutron fraction. And $\lambda_{i}$ is the decay constant of fission
product isotope $Q_{i}$.

The entire dynamic process of nuclear reactor is briefly described as follows.
It is not difficult to draw that if the initial value of neutrons is $N_{t}$,
neutron population through a nuclear fission is $kN_{t}$ where delayed neutron
density is $k\beta N_{t}$ and the number of new neutrons is $k(1-\beta)N_{t}$.
At first, discuss the change of $i$-th delayed neutron precursors density
$Q_{i}(t)$ with respect to time interval $\Delta t$. When a fission occurs,
this reactor produces delayed neutron precursors $Q_{i}$ with a certain rate
$\frac{k\beta_{i}N_{t}}{l}$ . Another situation, $i$-th delayed neutron
precursor decays into a neutron with decay constant $\lambda_{i}$. Then the
change of $Q_{i}$ is

$\Delta
Q_{i}=\bigg{(}-\lambda_{i}Q_{i}+\frac{k\beta_{i}N_{t}}{l}\bigg{)}\Delta t.$
(38)

Next, consider the change of the neutron population $N(t)$ with respect to
time interval $\Delta t$. There are three situations about the change of the
neutron population which are source, transformation and born. The source event
represents that an extraneous neutron source delivers $B$ neutrons to this
reactor at a constant rate, so derive $B\Delta t$ neutrons at time interval
$\Delta t$. The transformation event represents that $i$-th delayed neutron
precursor decays into a neutron with decay constant $\lambda_{i}$. The born
event represents neutrons produced when a fission occurs with rate
$\frac{k(1-\beta)-1}{l}N_{t}$. Given the above analysis, we can get the change
of $N(t)$:

$\Delta
N_{t}=\bigg{(}B+\sum_{i=1}^{6}\lambda_{i}Q_{i}+\frac{k(1-\beta)-1}{l}N_{t}\bigg{)}\Delta
t.$ (39)

According to Eq.(39) and Eq.(38), we can derive the following formula about
$N(t)$ and $Q_{i}(t)$

$\left\\{\begin{array}[]{l}\displaystyle\frac{\mathrm{d}N_{t}}{\mathrm{d}t}=B+\sum\limits_{i=1}^{6}\lambda_{i}Q_{i}+\frac{k(1-\beta)-1}{l}N_{t}\\\
\displaystyle\frac{\mathrm{d}Q_{i}}{\mathrm{d}t}=-\lambda_{i}Q_{i}+\frac{k\beta_{i}}{l}N_{t},\quad
i=1,2,\cdots,6.\end{array}\right.$ (40)

The six groups of delayed neutron dynamics Eqs. (40) are still complicated in
form. Considering a relatively stable situation, assume just one group of
delayed neutron precursor in this reactor. We can get $\lambda=\lambda_{1}$
and $\beta=\beta_{1}$. In this way, the seven equations are reduced to two
equations, so the equations are a little easier to solve.

In the beginning, the neutron population is very low in the reactor. In order
to accelerate the start-up speed, it is necessary to add an extraneous neutron
source B. As the reaction progresses, neutron population is much greater than
the source neutron concentration, so we can regard as $B=0$. Thus we can get
the equation of just one delayed neutron precursor:

$\left\\{\begin{array}[]{l}\vspace{5pt}\displaystyle\frac{\mathrm{d}N_{t}}{\mathrm{d}t}=\lambda
Q_{t}+\frac{k(1-\beta)-1}{l}N_{t}\\\
\displaystyle\frac{\mathrm{d}Q_{t}}{\mathrm{d}t}=-\lambda
Q_{t}+\frac{k\beta}{l}N_{t}.\end{array}\right.$ (41)

The nuclear reactor kinetics of one delayed neutron can roughly reflect the
fast transient of neutron population change and the slow transient with stable
period. Since we care more about the neutron population and the properties of
this system of Eq.(41), we first consider transforming it into a second-order
equation

$\frac{\mathrm{d}^{2}N_{t}}{\mathrm{d}t^{2}}=\bigg{(}\frac{k(1-\beta)-1}{l}-\lambda\bigg{)}\frac{\mathrm{d}N_{t}}{\mathrm{d}t}+\lambda\cdot\frac{k-1}{l}N_{t}.$
(42)

The actual expression of the above Eq.(42) is that the neutron population
varies with time in the same way. However, such a situation will only occur in
ideal conditions, or in computer simulations. Similar to the situation of
start-up, shutdown in reactor or discovery of unknown new radioactive
elements, we must fully consider the impact of unknown situations on the
reactor to ensure the safe operation of the reactor. Since in the actual
reaction, $\beta$ and $\lambda$ will be disturbed by environmental factors and
will not be a stable constant, that is

$\beta(t)=\beta+\sigma_{1}\cdot``\text{noise}"$
$\lambda(t)=\lambda+\sigma_{2}\cdot``\text{noise}".$

Based on uncertainty theory, the “noise” can be consider as a normal uncertain
variable $\mathcal{N}(0,1)$. In other words,

$\frac{\mathrm{d}C_{1t}}{\mathrm{d}t}=\frac{C_{1,t+\Delta t}-C_{1t}}{\Delta
t}$

and

$\frac{\mathrm{d}C_{2t}}{\mathrm{d}t}=\frac{C_{2,t+\Delta t}-C_{2t}}{\Delta
t}$

where $C_{1t}$ and $C_{2t}$ are Liu processes, so we get

$\left\\{\begin{array}[]{l}\displaystyle\beta(t)=\beta+\sigma_{1}\cdot\frac{\mathrm{d}C_{1t}}{\mathrm{d}t}\vspace{5pt}\\\
\displaystyle\lambda(t)=\lambda+\sigma_{2}\cdot\frac{\mathrm{d}C_{2t}}{\mathrm{d}t}\end{array}\right.$

where $\sigma_{1}$ and $\sigma_{2}$ are nonnegative constants, and they
represent the noise levels. Consequently, we derive the second-order uncertain
differential equation of nuclear reactor kinetics, i.e.

$\left\\{\begin{array}[]{l}\vspace{1ex}\displaystyle\frac{\mathrm{d}^{2}N_{t}}{\mathrm{d}t^{2}}\\!=\\!\Big{(}\frac{k(1\\!-\\!\beta)\\!-\\!1}{l}\\!-\\!\lambda\Big{)}\frac{\mathrm{d}N_{t}}{\mathrm{d}t}\\!+\\!\lambda\frac{k\\!-\\!\\!1}{l}N_{t}\\!-\\!\frac{k}{l}\sigma_{1}\frac{\mathrm{d}N_{t}}{\mathrm{d}t}\frac{\mathrm{d}C_{1t}}{\mathrm{d}t}\\\
\ \ \ \ \ \ \ \ \ \ \ \ \
\displaystyle+\sigma_{2}\big{(}\frac{k\\!-\\!1}{l}N_{t}\\!-\\!\frac{\mathrm{d}N_{t}}{\mathrm{d}t}\big{)}\frac{\mathrm{d}C_{2t}}{\mathrm{d}t}\\\
\vspace{1ex}N(t_{0})=n_{0}\\\
\displaystyle\frac{\mathrm{d}N_{t}}{\mathrm{d}t}\bigg{|}_{t=t_{0}}=n_{0}^{\prime}.\end{array}\right.$
(43)

###  6.2 Parameter estimation and hypothesis test

Although we deduced the uncertain nuclear reactor kinetics Eq.(43), only in
connection with the actual data can we reflect whether this equation is
reasonable or not. In this part, we will utilize experimental data to estimate
the unknown parameters in this equation. To safeguard real data, the utilized
data are based on specific parameters of an actual reactor.

Regard $\lambda=0.0785/sec,\beta=0.0065,k=1.001,l=10^{-4}sec$ as a numerical
example, which are the data involved in the thermal fission in a nuclear
reactor about uranium-$235$ fuel.

Consider the updated uncertain differential equation of nuclear reactor
kinetics

$\left\\{\begin{array}[]{l}\vspace{1ex}\displaystyle\frac{\mathrm{d}^{2}N_{t}}{\mathrm{d}t^{2}}\\!=\\!-55.1435\frac{\mathrm{d}N_{t}}{\mathrm{d}t}\\!+\\!0.785N_{t}\\!-\\!10010\sigma_{1}\frac{\mathrm{d}N_{t}}{\mathrm{d}t}\frac{\mathrm{d}C_{1t}}{\mathrm{d}t}\\\
\ \ \ \ \ \ \ \ \ \ \ \
\displaystyle+\sigma_{2}\bigg{(}10N_{t}\\!-\\!\frac{\mathrm{d}N_{t}}{\mathrm{d}t}\bigg{)}\frac{\mathrm{d}C_{2t}}{\mathrm{d}t}\\\
\vspace{1ex}N_{t_{j}}=x_{t_{j}}\vspace{1ex}\\\
\displaystyle\frac{\mathrm{d}N_{t}}{\mathrm{d}t}\bigg{|}_{t=t_{j}}=x^{\prime}_{t_{j}}\\\
\end{array}\right.$ (44)

where $x_{t_{j}}$, $x^{\prime}_{t_{j}}$ are observations at time $t_{j}$, and
$\sigma_{1}$ and $\sigma_{2}$ are unknown parameters to be estimated.

Table 3 shows $61$ observed data of this unclear reactor based on the values
given above. And the fluctuation of the neutron population over time is shown
in the Fig.5.

Figure 5: Neutron population from $t=0$ to $t=6$ with initial value
$n_{0}=1.2157$

Next, we will apply residuals and moment estimation to estimate the unknown
parameters. According to the Algorithm 1, we can produce 60 residuals

$\varepsilon_{1}(\sigma_{1},\sigma_{2}),\
\varepsilon_{2}(\sigma_{1},\sigma_{2}),\ \cdots,\
\varepsilon_{60}(\sigma_{1},\sigma_{2})$

for any given parameters $\sigma_{1}$ and $\sigma_{2}$.

Table 3: 61 experimental data of nuclear reactor

$t$ | $N_{t}$ | $t$ | $N_{t}$ | $t$ | $N_{t}$ | $t$ | $N_{t}$ | $t$ | $N_{t}$ | $t$ | $N_{t}$
---|---|---|---|---|---|---|---|---|---|---|---
$0$ | $1.2157$ | $1$ | $1.2313$ | $2$ | $1.2530$ | $3$ | $1.2658$ | $4$ | $1.2902$ | $5$ | $1.3249$
$0.1$ | $1.2165$ | $1.1$ | $1.2388$ | $2.1$ | $1.2545$ | $3.1$ | $1.2689$ | $4.1$ | $1.2922$ | $5.1$ | $1.3138$
$0.2$ | $1.2186$ | $1.2$ | $1.2339$ | $2.2$ | $1.2498$ | $3.2$ | $1.2717$ | $4.2$ | $1.2817$ | $5.2$ | $1.3064$
$0.3$ | $1.2213$ | $1.3$ | $1.2386$ | $2.3$ | $1.2572$ | $3.3$ | $1.2748$ | $4.3$ | $1.2849$ | $5.3$ | $1.3136$
$0.4$ | $1.2236$ | $1.4$ | $1.2361$ | $2.4$ | $1.2550$ | $3.4$ | $1.2797$ | $4.4$ | $1.2908$ | $5.4$ | $1.3047$
$0.5$ | $1.2236$ | $1.5$ | $1.2399$ | $2.5$ | $1.2656$ | $3.5$ | $1.2782$ | $4.5$ | $1.2964$ | $5.5$ | $1.3288$
$0.6$ | $1.2262$ | $1.6$ | $1.2420$ | $2.6$ | $1.2607$ | $3.6$ | $1.2880$ | $4.6$ | $1.3057$ | $5.6$ | $1.3241$
$0.7$ | $1.2271$ | $1.7$ | $1.2459$ | $2.7$ | $1.2697$ | $3.7$ | $1.2823$ | $4.7$ | $1.2975$ | $5.7$ | $1.3250$
$0.8$ | $1.2296$ | $1.8$ | $1.2492$ | $2.8$ | $1.2635$ | $3.8$ | $1.2922$ | $4.8$ | $1.3101$ | $5.8$ | $1.3145$
$0.9$ | $1.2317$ | $1.9$ | $1.2494$ | $2.9$ | $1.2688$ | $3.9$ | $1.2815$ | $4.9$ | $1.3050$ | $5.9$ | $1.3243$
|  |  |  |  |  |  |  |  |  | $6$ | $1.3285$

According to Eq.(34), the generalized moment estimation
($\hat{\sigma_{1}},\hat{\sigma_{2}}$) for ($\sigma_{1},\sigma_{2}$) is the
optimal solution of

$\left\\{\begin{array}[]{ll}\vspace{1mm}\displaystyle\min\limits_{\sigma_{1},\sigma_{2}}\sum\limits_{q=1}^{2}\bigg{(}\frac{1}{60}\sum_{i=1}^{60}\varepsilon_{i}^{q}(\sigma_{1},\sigma_{2})-\frac{1}{q+1}\bigg{)}^{2}\\\
\rm{subject\ to:}\\\ \
\qquad\qquad\sigma_{1},\sigma_{2}\in(0,1).\end{array}\right.$ (45)

Solving the above minimization problem (45) by MATLAB, we can obtain

$\hat{\sigma_{1}}=0.000143,\qquad\hat{\sigma_{2}}=0.296798,$

and the minimum value of Eq.(45) is

$\min_{\sigma_{1},\sigma_{2}}\sum_{q=1}^{2}\bigg{(}\frac{1}{60}\sum_{i=1}^{60}\varepsilon_{i}^{q}(\sigma_{1},\sigma_{2})-\frac{1}{q+1}\bigg{)}^{2}=8.25350\times
10^{-12}.$

Thus we obtain uncertain differential equation of this nuclear reactor

$\left\\{\begin{array}[]{ll}\vspace{2ex}\displaystyle\frac{\mathrm{d}^{2}N_{t}}{\mathrm{d}t^{2}}=\displaystyle-55.1435\frac{\mathrm{d}N_{t}}{\mathrm{d}t}\\!+\\!0.785N_{t}\\!-\\!1.43143\frac{\mathrm{d}N_{t}}{\mathrm{d}t}\cdot\frac{\mathrm{d}C_{1t}}{\mathrm{d}t}\\!\\\
\ \ \ \ \ \ \ \ \ \ \ \ \
\displaystyle+\bigg{(}2.96798N_{t}\\!-\\!0.296798\frac{\mathrm{d}N_{t}}{\mathrm{d}t}\bigg{)}\frac{\mathrm{d}C_{2t}}{\mathrm{d}t}\\\
\vspace{2ex}N(0)=1.2157\\\ N^{\prime}(0)=0.008.\end{array}\right.$ (46)

Moreover, we can also get 60 residuals
$\varepsilon_{1},\cdots,\varepsilon_{60}$ of uncertain nuclear reactor
kinetics Eq.(46) as shown in Table 4. On the basis of hypothesis test, we need
to test whether these 60 residuals $\varepsilon_{1},\cdots,\varepsilon_{60}$
have a good fit to the linear uncertainty distribution $\mathcal{L}(0,1)$. The
test method is as follows.

Table 4: 60 residuals

$j$ | $\varepsilon_{j}$ | $j$ | $\varepsilon_{j}$ | $j$ | $\varepsilon_{j}$ | $j$ | $\varepsilon_{j}$ | $j$ | $\varepsilon_{j}$ | $j$ | $\varepsilon_{j}$
---|---|---|---|---|---|---|---|---|---|---|---
$1$ | $0.4373$ | $11$ | $0.8227$ | $21$ | $0.4787$ | $31$ | $0.5851$ | $41$ | $0.5070$ | $51$ | $0.0317$
$2$ | $0.5238$ | $12$ | $0.1370$ | $22$ | $0.1454$ | $32$ | $0.5664$ | $42$ | $0.0342$ | $52$ | $0.0809$
$3$ | $0.5632$ | $13$ | $0.6916$ | $23$ | $0.8153$ | $33$ | $0.5845$ | $43$ | $0.5906$ | $53$ | $0.7937$
$4$ | $0.5436$ | $14$ | $0.2351$ | $24$ | $0.2537$ | $34$ | $0.6861$ | $44$ | $0.7404$ | $54$ | $0.0540$
$5$ | $0.3800$ | $15$ | $0.6358$ | $25$ | $0.9065$ | $35$ | $0.2970$ | $45$ | $0.7223$ | $55$ | $0.9952$
$6$ | $0.5616$ | $16$ | $0.5177$ | $26$ | $0.1407$ | $36$ | $0.8839$ | $46$ | $0.8696$ | $56$ | $0.1549$
$7$ | $0.4387$ | $17$ | $0.6417$ | $27$ | $0.8673$ | $37$ | $0.1187$ | $47$ | $0.0638$ | $57$ | $0.4357$
$8$ | $0.5498$ | $18$ | $0.6019$ | $28$ | $0.1011$ | $38$ | $0.8856$ | $48$ | $0.9360$ | $58$ | $0.0368$
$9$ | $0.5226$ | $19$ | $0.3940$ | $29$ | $0.7186$ | $39$ | $0.0329$ | $49$ | $0.1401$ | $59$ | $0.8784$
$10$ | $0.3566$ | $20$ | $0.6192$ | $30$ | $0.2134$ | $40$ | $0.8542$ | $50$ | $0.9875$ | $60$ | $0.6384$

Given a significance level $\alpha=0.05$, and due to $\alpha\times 60=3$, the
test is

$\displaystyle W=\\{(\varepsilon_{1},\cdots,\varepsilon_{60}):$
$\displaystyle{\text{there are at least 3 of index j's with }}1\leq j\leq 60$
$\displaystyle{\text{ such that }}\varepsilon_{j}<0.025{\text{ or
}}\varepsilon_{j}>0.975\\}.$

From Table 4 and Fig.6, it is clearly that $\varepsilon_{50}$,
$\varepsilon_{55}$ are the only two residuals not in $[0.025,0.975]$.
Therefore, $(\varepsilon_{1},\cdots,\varepsilon_{60})\notin W$. In other
words, Eq.(46) really have a good fit with the observed data.

Figure 6: Residual plot of hypothesis test

Why do we use uncertain differential equations to describe nuclear reactor
kinetics? The reason is as follows. When we divide the 60 residuals into

$(\varepsilon_{1},\varepsilon_{2},\cdots,\varepsilon_{21})\ \mbox{and}\
(\varepsilon_{44},\varepsilon_{45},\cdots,\varepsilon_{60}),$

the two-sample Kolmogorov-Smirnov test showed that the above two parts from
the residuals do not come from the same population via the function ”kstest2”
in Matlab. Thus the residuals $\varepsilon_{1},\cdots,\varepsilon_{60}$ are
not white noise in the sense of probability theory.

Next, consider the $\alpha$-path of the uncertain nuclear reactor kinetics
Eq.(46). $N_{t}^{\alpha}$ is the solution of

$\left\\{\begin{array}[]{l}\displaystyle\frac{\mathrm{d}^{2}N_{t}^{\alpha}}{\mathrm{d}t^{2}}=\displaystyle-55.1435\frac{\mathrm{d}N_{t}^{\alpha}}{\mathrm{d}t}\\!+\\!0.785N_{t}^{\alpha}\\!+\bigg{|}-\\!1.43143\frac{\mathrm{d}N_{t}^{\alpha}}{\mathrm{d}t}\bigg{|}\Phi^{-1}(\alpha)\\\
\ \ \ \ \ \ \ \ \ \ \ \ \ \
\displaystyle+\bigg{|}2.96798N_{t}^{\alpha}\\!-\\!0.296798\frac{\mathrm{d}N_{t}^{\alpha}}{\mathrm{d}t}\bigg{|}\Phi^{-1}(\alpha)\\\
\vspace{2ex}N(0)=1.2157\\\ N^{\prime}(0)=0.008.\end{array}\right.$ (47)

In an actual controlled reactor, as the nuclear reaction goes on, the value of
$\frac{\mathrm{d}N_{t}}{\mathrm{d}t}/N_{t}$ is very small, so that $N_{t}$ is
much bigger than $\frac{\mathrm{d}N_{t}}{\mathrm{d}t}$. We can transform
Eq.(47) into

$\left\\{\begin{array}[]{l}\vspace{2ex}\displaystyle\frac{\mathrm{d}^{2}N_{t}^{\alpha}}{\mathrm{d}t^{2}}=\displaystyle\Big{(}-55.1435+1.134632\Phi^{-1}(\alpha)\Big{)}\frac{\mathrm{d}N_{t}^{\alpha}}{\mathrm{d}t}\\\
\ \ \ \ \ \ \ \ \ \ \ \ \
\displaystyle+\Big{(}0.785+2.96798\Phi^{-1}(\alpha)\Big{)}N_{t}^{\alpha}\\\
\vspace{2ex}N(0)=1.2157\\\ N^{\prime}(0)=0.008.\end{array}\right.$ (48)

When $\alpha>0.4$,

$\Big{(}-55.1435+1.134632\Phi^{-1}(\alpha)\Big{)}\frac{\mathrm{d}N_{t}^{\alpha}}{\mathrm{d}t}+\Big{(}0.785+2.96798\Phi^{-1}(\alpha)\Big{)}N_{t}^{\alpha}$

is a monotonically increasing function with respect to $N_{t}^{\alpha}$. Thus,
by Theorem 2.3, for the uncertain nuclear reactor kinetics Eq.(46), the
inverse uncertainty distribution $\Psi_{t}^{-1}(\alpha)$ is the solution
$N_{t}^{\alpha}$ of Eq.(48), i.e.,

$\begin{array}[]{l}\displaystyle\Psi_{t}^{-1}(\alpha)\\!=\\!\mathrm{exp}\left\\{\frac{(P\\!-\\!\sqrt{T})}{2}t\right\\}\cdot\frac{-0.016\\!+\\!1.2157(\sqrt{T}\\!+\\!P)}{2\sqrt{T}}\\\
\ \ \ \ \ \ \ \ \ \ \ \ \ \
\displaystyle+\mathrm{exp}\left\\{\frac{(P\\!+\\!\sqrt{T})}{2}t\right\\}\cdot\frac{0.016\\!+\\!1.2157(\sqrt{T}\\!-\\!P)}{2\sqrt{T}}\end{array}$

where

$\displaystyle P=-55.1435+1.134632\Phi^{-1}(\alpha),$ $\displaystyle
R=0.785+2.96798\Phi^{-1}(\alpha),$ $\displaystyle T=P^{2}+4R,$
$\displaystyle\Phi^{-1}(\alpha)=\frac{\sqrt{3}}{\pi}\mathrm{ln}\frac{\alpha}{1-\alpha},\quad
0.4<\alpha<1.$

##  7 Conclusion

The comparison theorems of high-order ordinary differential equations are
rigorously proved. And the $\alpha$-path of high-order uncertain differential
equation is presented. The method to solve a family of high-order uncertain
differential equations is proposed including parameter estimation and
hypothesis test. Uncertain nuclear reactor kinetics equation is introduced as
an example to illustrate this method.

## References

  * [1] Allen, E.(1999). Stochastic differential equations and persistence time for two interacting populations. Dynamics of Continuous Discrete and Impulsive Systems, 5(1-4), 271-281.
  * [2] Hayes, J. G., & Allen, E. J. (2005). Stochastic point-kinetics equations in nuclear reactor dynamics. Annals of Nuclear Energy, 32(6), 572-587.
  * [3] Kinard, M., & Allen, E. J.(2004). Efficient numerical solution of the point kinetics equations in nuclear reactor dynamics. Annals of Nuclear Energy, 31(9), 1039-1051.
  * [4] Lio, W., & Liu, B.(2021). Initial value estimation of uncertain differential equations and zero-day of COVID-19 spread in China. Fuzzy Optimization and Decision Making, 20(2), 177–188.
  * [5] Liu, B. (2007). Uncertainty Theory (2nd ed.). Berlin, Germany: Springer-Verlag.
  * [6] Liu, B. (2008). Fuzzy process, hybrid process and uncertain process. Journal of Uncertain Systems, 2(1), 3-16.
  * [7] Liu, B. (2015). Uncertainty Theory (4th ed.). Berlin: Springer.
  * [8] Liu, Y., & Liu, B. (2022). Residual analysis and parameter estimation of uncertain differential equations. Fuzzy Optimization and Decision Making, 21(4), 513-530.
  * [9] Liu, Y., & Liu, B. (2023a). A modified uncertain maximum likelihood estimation with applications in uncertain statistics. Communications in Statistics - Theory and Methods. https://doi.org/10.1080/03610926.2023.2248534.
  * [10] Liu, Y., & Liu, B. (2023b). Estimation of uncertainty distribution function by the principle of least squares. Communications in Statistics - Theory and Methods. https://doi.org/10.1080/03610926.2023.2269451.
  * [11] Liu, Z., & Yang, Y. (2021). Uncertain pharmacokinetic model based on uncertain differential equation. Applied Mathematics and Computation, 404, 126118.
  * [12] Mehrdoust, F., Noorani, I., & Xu, W. (2023). Uncertain energy model for electricity and gas futures with application in spark-spread option price. Fuzzy Optimization and Decision Making, 22(1), 123-148.
  * [13] Tang, H., & Yang, X. (2021). Uncertain chemical reaction equation. Applied Mathematics and Computation, 411, 126479.
  * [14] Yang, L., & Liu, Y. (2023). Solution method and parameter estimation of uncertain partial differential equation with application to China’s population. Fuzzy Optimization and Decision Making, 23(1), 155-177.
  * [15] Yang, X., & Ke, H. (2023). Uncertain interest rate model for Shanghai interbank offered rate and pricing of American swaption. Fuzzy Optimization and Decision Making, 22(3), 447–462.
  * [16] Yao, K. (2016). Uncertainty differential equations. Berlin: Springer-Verlag.
  * [17] Yao, K., & Chen, X. (2013). A numerical method for solving uncertain differential equations. Journal of Intelligent and Fuzzy Systems, 25(3), 825-832.
  * [18] Ye, T. (2023). Partial derivatives of uncertain fields and uncertain partial differential equations. Fuzzy Optimization and Decision Making, https://doi.org/10.1007/s10700-023-09417-3.
  * [19] Ye, T., & Liu, B. (2023). Uncertain hypothesis test for uncertain differential equations. Fuzzy Optimization and Decision Making, 22(2), 195-211.
  * [20] Ye, T., & Zheng, H. (2023). Analysis of birth rates in China with uncertain statistics. Journal of intelligent and fuzzy systems, 44(6), 10621-10632.
  * [21] Zhu, Y. (2023). On uncertain partial differential equations. Fuzzy Optimization and Decision Making, https://doi.org/10.1007/s10700-023-09418-2.


Table 8: Momentum Lookback Window This table shows the annualized out-of-sample Sharpe Ratio of various portfolio selection methods based on momentum signals of different lookback window sizes. Panels A-E display the results for 42 days, 63 days, 84 days, 126 days, and 252 days. The test assets are the 25 and 100 equity portfolios independently sorted on size and three other factors (book-to-market, operating profitability, and investment). The numbers in parenthesis are the t-statistics for the portfolio’s Sharpe ratio and information ratio using the standard error of Lo, (2002). The covariances are estimated on a rolling walk-forward basis with $T=120$ non-overlapping monthly returns. All measures are computed with 579 monthly portfolio weights from 09-09-1974 until 10-21-2022.  | FF25 | FF100 |  ME/OP 25 |  ME/OP 100 |  ME/INV 25 |  ME/INV 100
---|---|---|---|---|---|---
Panel A. Sharpe Ratio, Lookback Window Size: 42
CP2 | $\underset{(6.054)}{0.887}$ | $\underset{(6.837)}{1.006}$ | $\underset{(5.895)}{0.863}$ | $\underset{(5.074)}{0.740}$ | $\underset{(6.198)}{0.909}$ | $\underset{(4.853)}{0.707}$
PP2 | $\underset{(4.465)}{0.649}$ | $\underset{(4.298)}{0.624}$ | $\underset{(4.587)}{0.668}$ | $\underset{(3.759)}{0.545}$ | $\underset{(5.342)}{0.780}$ | $\underset{(4.728)}{0.688}$
MVO | $\underset{(1.647)}{0.238}$ | $\underset{(1.380)}{0.199}$ | $\underset{(1.444)}{0.208}$ | $\underset{(-0.355)}{-0.051}$ | $\underset{(1.468)}{0.212}$ | $\underset{(0.761)}{0.110}$
UNI | $\underset{(3.031)}{0.438}$ | $\underset{(2.532)}{0.366}$ | $\underset{(3.462)}{0.502}$ | $\underset{(2.748)}{0.398}$ | $\underset{(2.982)}{0.432}$ | $\underset{(2.265)}{0.327}$
Panel B. Sharpe Ratio, Lookback Window Size: 63
CP2 | $\underset{(4.575)}{0.665}$ | $\underset{(6.681)}{0.982}$ | $\underset{(5.948)}{0.871}$ | $\underset{(5.389)}{0.787}$ | $\underset{(6.143)}{0.901}$ | $\underset{(6.199)}{0.909}$
PP2 | $\underset{(4.216)}{0.612}$ | $\underset{(3.962)}{0.575}$ | $\underset{(4.913)}{0.716}$ | $\underset{(4.555)}{0.663}$ | $\underset{(4.627)}{0.674}$ | $\underset{(4.479)}{0.652}$
MVO | $\underset{(2.273)}{0.328}$ | $\underset{(1.995)}{0.288}$ | $\underset{(1.284)}{0.185}$ | $\underset{(-1.354)}{-0.195}$ | $\underset{(1.753)}{0.253}$ | $\underset{(1.229)}{0.177}$
UNI | $\underset{(2.901)}{0.420}$ | $\underset{(2.536)}{0.366}$ | $\underset{(3.308)}{0.479}$ | $\underset{(2.325)}{0.336}$ | $\underset{(2.716)}{0.393}$ | $\underset{(1.979)}{0.286}$
Panel C. Sharpe Ratio, Lookback Window Size: 84
CP2 | $\underset{(4.306)}{0.625}$ | $\underset{(6.196)}{0.908}$ | $\underset{(4.628)}{0.674}$ | $\underset{(3.743)}{0.543}$ | $\underset{(4.237)}{0.616}$ | $\underset{(5.834)}{0.854}$
PP2 | $\underset{(3.504)}{0.508}$ | $\underset{(3.129)}{0.453}$ | $\underset{(2.826)}{0.409}$ | $\underset{(3.329)}{0.482}$ | $\underset{(3.211)}{0.465}$ | $\underset{(3.460)}{0.502}$
MVO | $\underset{(1.829)}{0.264}$ | $\underset{(2.110)}{0.305}$ | $\underset{(1.117)}{0.161}$ | $\underset{(-0.675)}{-0.097}$ | $\underset{(1.480)}{0.214}$ | $\underset{(2.169)}{0.313}$
UNI | $\underset{(1.917)}{0.277}$ | $\underset{(1.891)}{0.273}$ | $\underset{(2.352)}{0.340}$ | $\underset{(1.796)}{0.259}$ | $\underset{(1.715)}{0.248}$ | $\underset{(1.678)}{0.242}$
Panel D. Sharpe Ratio, Lookback Window Size: 126
CP2 | $\underset{(4.605)}{0.670}$ | $\underset{(6.042)}{0.885}$ | $\underset{(3.755)}{0.545}$ | $\underset{(3.352)}{0.486}$ | $\underset{(4.179)}{0.607}$ | $\underset{(5.270)}{0.769}$
PP2 | $\underset{(3.169)}{0.459}$ | $\underset{(2.957)}{0.428}$ | $\underset{(3.055)}{0.442}$ | $\underset{(3.865)}{0.561}$ | $\underset{(2.486)}{0.359}$ | $\underset{(1.907)}{0.275}$
MVO | $\underset{(2.172)}{0.314}$ | $\underset{(2.620)}{0.379}$ | $\underset{(0.808)}{0.117}$ | $\underset{(-0.743)}{-0.107}$ | $\underset{(2.331)}{0.337}$ | $\underset{(1.780)}{0.257}$
UNI | $\underset{(2.455)}{0.355}$ | $\underset{(2.154)}{0.311}$ | $\underset{(2.390)}{0.346}$ | $\underset{(1.910)}{0.276}$ | $\underset{(2.040)}{0.295}$ | $\underset{(1.760)}{0.254}$
Panel E. Sharpe Ratio, Lookback Window Size: 252
CP2 | $\underset{(6.851)}{1.008}$ | $\underset{(8.089)}{1.200}$ | $\underset{(5.511)}{0.805}$ | $\underset{(4.953)}{0.722}$ | $\underset{(6.460)}{0.949}$ | $\underset{(7.629)}{1.129}$
PP2 | $\underset{(4.035)}{0.585}$ | $\underset{(3.526)}{0.511}$ | $\underset{(3.683)}{0.534}$ | $\underset{(3.654)}{0.530}$ | $\underset{(3.521)}{0.511}$ | $\underset{(2.282)}{0.330}$
MVO | $\underset{(3.458)}{0.501}$ | $\underset{(3.986)}{0.578}$ | $\underset{(2.502)}{0.362}$ | $\underset{(0.300)}{0.043}$ | $\underset{(3.240)}{0.469}$ | $\underset{(2.701)}{0.391}$
UNI | $\underset{(3.023)}{0.437}$ | $\underset{(2.672)}{0.386}$ | $\underset{(3.612)}{0.524}$ | $\underset{(3.103)}{0.449}$ | $\underset{(3.042)}{0.441}$ | $\underset{(2.650)}{0.383}$

Table 9: Approximate versus Actual Solution This table shows the out-of-sample performances (in percentages) of the ‘approximated’ versus ‘actual’ portfolio policy using a 21-day momentum signal. Panels A-G display the results for the average returns, the Sharpe ratio of returns, the alpha with respect to the Fama-French 5-factor regression augmented with the univariate factor (UNI), the beta of the simple factor, the idiosyncratic volatility, and the information ratio. The test assets are the 25 and 100 equity portfolios independently sorted on size and three other factors (book-to-market, operating profitability, and investment). The numbers in parenthesis are the t-statistics for the portfolio’s Sharpe ratio and information ratio using the standard error of Lo, (2002). The covariances are estimated on a rolling walk-forward basis with $T=120$ non-overlapping monthly returns. All measures are computed with 579 monthly portfolio weights from 09-09-1974 until 10-21-2022. | FF25 | FF100 |  ME/OP 25 |  ME/OP 100 |  ME/INV 25 |  ME/INV 100
---|---|---|---|---|---|---
Panel A. Average Returns (%)
CP2 (Approx) | $2.199$ | $1.341$ | $1.731$ | $0.957$ | $1.797$ | $0.777$
CP2 (Full) | $2.253$ | $1.175$ | $1.578$ | $0.822$ | $1.675$ | $0.596$
Panel B. Standard Deviation of Returns (%)
CP2 (Approx) | $2.177$ | $1.366$ | $2.049$ | $1.260$ | $2.057$ | $1.248$
CP2 (Full) | $2.210$ | $1.356$ | $2.119$ | $1.229$ | $2.061$ | $1.309$
Panel C. Sharpe Ratio
CP2 (Approx) | $\underset{(6.872)}{1.011}$ | $\underset{(6.685)}{0.982}$ | $\underset{(5.797)}{0.848}$ | $\underset{(5.226)}{0.763}$ | $\underset{(5.985)}{0.877}$ | $\underset{(4.308)}{0.626}$
CP2 (Full) | $\underset{(6.933)}{1.020}$ | $\underset{(5.934)}{0.868}$ | $\underset{(5.140)}{0.750}$ | $\underset{(4.620)}{0.672}$ | $\underset{(5.588)}{0.817}$ | $\underset{(3.177)}{0.460}$
Panel D. Alpha (%)
CP2 (Approx) | $1.849$ | $1.197$ | $1.109$ | $0.826$ | $1.334$ | $0.679$
CP2 (Full) | $1.908$ | $1.043$ | $0.942$ | $0.739$ | $1.224$ | $0.536$
Panel E. Beta to Simple Factor
CP2 (Approx) | $0.146$ | $0.071$ | $0.259$ | $0.086$ | $0.172$ | $0.053$
CP2 (Full) | $0.144$ | $0.067$ | $0.267$ | $0.057$ | $0.165$ | $0.046$
Panel F. Idiosyncratic Volatility (%)
CP2 (Approx) | $2.043$ | $1.321$ | $1.749$ | $1.220$ | $1.911$ | $1.229$
CP2 (Full) | $2.083$ | $1.315$ | $1.807$ | $1.208$ | $1.921$ | $1.297$
Panel G. Information Ratio
CP2 (Approx) | $\underset{(6.176)}{0.905}$ | $\underset{(6.183)}{0.906}$ | $\underset{(4.360)}{0.634}$ | $\underset{(4.651)}{0.677}$ | $\underset{(4.792)}{0.698}$ | $\underset{(3.809)}{0.553}$
CP2 (Full) | $\underset{(6.247)}{0.916}$ | $\underset{(5.434)}{0.793}$ | $\underset{(3.596)}{0.522}$ | $\underset{(4.210)}{0.612}$ | $\underset{(4.379)}{0.637}$ | $\underset{(2.858)}{0.414}$

Table 10: Performance Summary for Two Signals Case This table shows the out-of-sample performances (in percentages) of various portfolio selection methods using two momentum signals with a lookback window of 21 days and 252 days. Panels A-G display the results for the average returns, the Sharpe ratio of returns, the alpha with respect to the Fama-French 5-factor regression augmented with the univariate factor (UNI), the beta of the simple factor, the idiosyncratic volatility, and the information ratio. The test assets are 25 and 100 equity portfolios independently sorted on size and three other factors (book-to-market, operating profitability, and investment). The numbers in parenthesis are the t-statistics for the portfolio’s Sharpe ratio and information ratio using the standard error of Lo, (2002). The covariances are estimated on a rolling walk-forward basis with $T=120$ non-overlapping monthly returns. All measures are computed with 579 monthly portfolio weights from 09-09-1974 until 10-21-2022. | FF25 | FF100 |  ME/OP 25 |  ME/OP 100 |  ME/INV 25 |  ME/INV 100
---|---|---|---|---|---|---
Panel A. Average Returns (%)
CP2 | $2.925$ | $1.807$ | $2.308$ | $1.027$ | $2.274$ | $1.551$
PP2 | $4.160$ | $4.143$ | $3.779$ | $4.082$ | $2.893$ | $2.522$
MVO | $0.761$ | $0.581$ | $0.595$ | $-0.032$ | $0.746$ | $0.386$
UNI | $2.882$ | $2.316$ | $2.743$ | $2.047$ | $2.431$ | $1.838$
Panel B. Standard Deviation of Returns (%)
CP2 | $2.343$ | $1.384$ | $2.153$ | $1.312$ | $2.072$ | $1.293$
PP2 | $6.042$ | $6.053$ | $5.288$ | $5.488$ | $4.707$ | $4.760$
MVO | $1.893$ | $1.431$ | $1.807$ | $1.332$ | $1.642$ | $1.347$
UNI | $5.489$ | $5.097$ | $4.543$ | $4.095$ | $4.619$ | $4.280$
Panel C. Sharpe Ratio
CP2 | $\underset{(8.365)}{1.244}$ | $\underset{(8.735)}{1.302}$ | $\underset{(7.266)}{1.072}$ | $\underset{(5.387)}{0.786}$ | $\underset{(7.426)}{1.096}$ | $\underset{(8.073)}{1.197}$
PP2 | $\underset{(4.849)}{0.706}$ | $\underset{(4.823)}{0.702}$ | $\underset{(5.005)}{0.729}$ | $\underset{(5.197)}{0.758}$ | $\underset{(4.338)}{0.630}$ | $\underset{(3.777)}{0.548}$
MVO | $\underset{(2.836)}{0.410}$ | $\underset{(2.851)}{0.412}$ | $\underset{(2.336)}{0.337}$ | $\underset{(-0.120)}{-0.017}$ | $\underset{(3.184)}{0.461}$ | $\underset{(2.028)}{0.293}$
UNI | $\underset{(3.766)}{0.546}$ | $\underset{(3.285)}{0.476}$ | $\underset{(4.263)}{0.619}$ | $\underset{(3.561)}{0.516}$ | $\underset{(3.751)}{0.544}$ | $\underset{(3.093)}{0.448}$
Panel D. Alpha (%)
CP2 | $2.287$ | $1.576$ | $1.508$ | $0.781$ | $1.641$ | $1.300$
PP2 | $2.048$ | $2.097$ | $1.691$ | $2.263$ | $1.318$ | $1.205$
MVO | $0.025$ | $0.185$ | $-0.195$ | $-0.391$ | $0.117$ | $-0.003$
Panel E. Beta to Simple Factor
CP2 | $0.206$ | $0.100$ | $0.277$ | $0.106$ | $0.202$ | $0.115$
PP2 | $0.735$ | $0.866$ | $0.783$ | $0.922$ | $0.600$ | $0.716$
MVO | $0.254$ | $0.174$ | $0.291$ | $0.171$ | $0.257$ | $0.198$
Panel F. Idiosyncratic Volatility (%)
CP2 | $2.040$ | $1.285$ | $1.745$ | $1.240$ | $1.856$ | $1.198$
PP2 | $4.494$ | $4.155$ | $3.909$ | $3.986$ | $3.858$ | $3.673$
MVO | $1.288$ | $1.124$ | $1.235$ | $1.142$ | $1.136$ | $1.058$
Panel G. Information Ratio
CP2 | $\underset{(7.584)}{1.121}$ | $\underset{(8.254)}{1.226}$ | $\underset{(5.906)}{0.864}$ | $\underset{(4.336)}{0.630}$ | $\underset{(6.037)}{0.884}$ | $\underset{(7.349)}{1.084}$
PP2 | $\underset{(3.149)}{0.456}$ | $\underset{(3.485)}{0.505}$ | $\underset{(2.990)}{0.433}$ | $\underset{(3.914)}{0.568}$ | $\underset{(2.364)}{0.342}$ | $\underset{(2.272)}{0.328}$
MVO | $\underset{(0.133)}{0.019}$ | $\underset{(1.139)}{0.164}$ | $\underset{(-1.095)}{-0.158}$ | $\underset{(-2.373)}{-0.343}$ | $\underset{(0.715)}{0.103}$ | $\underset{(-0.019)}{-0.003}$

Table 11: Weight Statistics for Two Signals Case This table shows the average monthly weight statistics of various portfolio selection methods estimated using a 21-day momentum signal. Panels A-F display the results for the turnover, the proportional leverage, the sum of negative weights in the portfolio, and the minimum and maximum portfolio weights. The test assets are the 25 and 100 equity portfolios independently sorted on size and three other factors (book-to-market, operating profitability, and investment). The covariances are estimated on a rolling walk-forward basis with $T=120$ non-overlapping monthly returns. All measures are computed with 579 monthly portfolio weights from 09-09-1974 until 10-21-2022.  | FF25 | FF100 |  ME/OP 25 |  ME/OP 100 |  ME/INV 25 |  ME/INV 100
---|---|---|---|---|---|---
Panel A. Turnover
CP2 | $0.749$ | $0.437$ | $0.806$ | $0.542$ | $0.762$ | $0.466$
PP2 | $0.806$ | $0.763$ | $0.878$ | $0.794$ | $0.891$ | $0.856$
MVO | $1.064$ | $1.075$ | $1.065$ | $1.082$ | $1.060$ | $1.079$
UNI | $0.959$ | $0.978$ | $0.970$ | $0.990$ | $0.969$ | $0.989$
Panel B. Proportional Leverage
CP2 | $0.477$ | $0.473$ | $0.480$ | $0.488$ | $0.470$ | $0.486$
PP2 | $0.478$ | $0.483$ | $0.472$ | $0.473$ | $0.477$ | $0.484$
MVO | $0.506$ | $0.499$ | $0.509$ | $0.503$ | $0.497$ | $0.497$
UNI | $0.485$ | $0.498$ | $0.488$ | $0.497$ | $0.485$ | $0.498$
Panel C. Sum of Negative Weights
CP2 | $-0.483$ | $-0.480$ | $-0.486$ | $-0.482$ | $-0.481$ | $-0.479$
PP2 | $-0.500$ | $-0.500$ | $-0.500$ | $-0.500$ | $-0.500$ | $-0.500$
MVO | $-0.499$ | $-0.499$ | $-0.500$ | $-0.501$ | $-0.498$ | $-0.498$
UNI | $-0.500$ | $-0.500$ | $-0.500$ | $-0.500$ | $-0.500$ | $-0.500$
Panel D. Minimum Weight
CP2 | $-0.121$ | $-0.035$ | $-0.113$ | $-0.036$ | $-0.114$ | $-0.035$
PP2 | $-0.104$ | $-0.034$ | $-0.108$ | $-0.037$ | $-0.104$ | $-0.033$
MVO | $-0.099$ | $-0.031$ | $-0.096$ | $-0.030$ | $-0.097$ | $-0.031$
UNI | $-0.090$ | $-0.025$ | $-0.089$ | $-0.025$ | $-0.089$ | $-0.024$
Panel E. Maximum Weight
CP2 | $0.109$ | $0.032$ | $0.102$ | $0.032$ | $0.102$ | $0.032$
PP2 | $0.092$ | $0.031$ | $0.096$ | $0.033$ | $0.096$ | $0.030$
MVO | $0.102$ | $0.031$ | $0.099$ | $0.031$ | $0.097$ | $0.030$
UNI | $0.088$ | $0.025$ | $0.087$ | $0.025$ | $0.089$ | $0.025$
Panel F. Standard Deviation of Weights
CP2 | $0.053$ | $0.013$ | $0.052$ | $0.013$ | $0.052$ | $0.013$
PP2 | $0.051$ | $0.013$ | $0.052$ | $0.013$ | $0.051$ | $0.013$
MVO | $0.051$ | $0.013$ | $0.050$ | $0.012$ | $0.050$ | $0.012$
UNI | $0.049$ | $0.012$ | $0.049$ | $0.012$ | $0.049$ | $0.012$

##  Appendix B  Proofs for all the Propositions

  * Proof of Proposition 1.

We provide the proof of this proposition in two parts.

Part 1. We drop the time subscript $t$ for brevity. From the cyclic property
of the trace operator and linearity of the expectation, the expected value of
the portfolio returns is

$\displaystyle\mathbb{E}[x^{\prime}Ar]=\mathbb{E}[\mathsf{Tr}(Arx^{\prime})]=\mathsf{Tr}(A\Sigma_{rx}).$
(42)

In order to analyze the variance of the portfolio returns, we appeal to the
following one-dimensional from Isserlis, (1918) or Wick, (1950), which
expresses the higher moments of centered multivariate Gaussian vectors in
terms of its second-moments.

######  Theorem 1.

Let $z_{1},z_{2},z_{3}$, and $z_{4}$ be jointly Gaussian random variables with
mean zero. Then we have the following results:

$\displaystyle\mathbb{E}[z_{1}]$ $\displaystyle=0\hskip 70.0pt$ (43)
$\displaystyle\mathbb{E}[z_{1}z_{2}]$
$\displaystyle=\mathsf{Cov}(z_{1},z_{2})$
$\displaystyle\mathbb{E}[z_{1}z_{2}z_{3}]$ $\displaystyle=0\hskip 70.0pt$
$\displaystyle\mathbb{E}[z_{1}z_{2}z_{3}z_{4}]$
$\displaystyle=\mathbb{E}[z_{1}z_{2}]\mathbb{E}[z_{3}z_{4}]+\mathbb{E}[z_{1}z_{3}]\mathbb{E}[z_{2}z_{4}]+\mathbb{E}[z_{1}z_{4}]\mathbb{E}[z_{2}z_{3}]$

The last equation of (43) takes all partitions of size two of the four
variables, which gives us three separate terms. By recasting fourth-order
terms using the covariance terms, we can express the variance as

$\displaystyle\mathsf{Var}[x^{\prime}Ar]$
$\displaystyle=\mathbb{E}[(x^{\prime}Ar)^{2}]-\mathbb{E}[x^{\prime}Ar]^{2}$
(44)
$\displaystyle=\sum_{i,j,k,l}\mathbb{E}[A_{ij}A_{kl}x_{i}r_{j}x_{k}r_{l}]-\sum_{i,j,k,l}\mathbb{E}[A_{ij}x_{i}r_{j}]\mathbb{E}[A_{kl}x_{k}r_{l}]$
$\displaystyle=\sum_{i,j,k,l}A_{ij}A_{kl}\mathbb{E}[x_{i}r_{j}]\mathbb{E}[x_{k}r_{l}]+\sum_{i,j,k,l}A_{ij}A_{kl}\mathbb{E}[x_{i}x_{k}]\mathbb{E}[r_{j}r_{l}]$
$\displaystyle\qquad+\sum_{i,j,k,l}A_{ij}A_{kl}\mathbb{E}[x_{i}r_{l}]\mathbb{E}[x_{k}r_{j}]-\sum_{i,j,k,l}A_{ij}A_{kl}\mathbb{E}[x_{i}r_{j}]\mathbb{E}[x_{k}r_{l}]$
$\displaystyle=\sum_{i,j,k,l}A_{ij}A_{kl}\mathbb{E}[x_{i}x_{k}]\mathbb{E}[r_{j}r_{l}]+\sum_{i,j,k,l}A_{ij}A_{kl}\mathbb{E}[x_{i}r_{l}]\mathbb{E}[x_{k}r_{j}]$
$\displaystyle=\sum_{i,j,k,l}\mathbb{E}[x_{k}x_{i}]A_{ij}\mathbb{E}[r_{j}r_{l}]A_{kl}+\sum_{i,j,k,l}\mathbb{E}[r_{l}x_{i}]A_{ij}\mathbb{E}[r_{j}x_{k}]A_{kl}.$

where $A_{ij}$ is the $(i,j)$-th entry of $A$. Reverting back to matrix
notation, we have

$\displaystyle\mathsf{Var}[x^{\prime}Ar]$
$\displaystyle=\mathsf{Tr}(\Sigma_{x}A\Sigma_{r}A^{\prime})+\mathsf{Tr}(\Sigma_{rx}A\Sigma_{rx}A).$
(45)

Part 2. In order to solve the full mean-variance objective function (4), we
use the tools from CCA that we have laid the ground in Section 3.1. To assist
in our derivation, we expand our notation and let $Q_{r}$ and $Q_{x}$ be
matrices whose columns contain the so-called $i$th canonical directions
$q_{r,i}:=\Sigma_{r}^{-1/2}u_{i}$ and $q_{x,i}:=\Sigma_{x}^{-1/2}v_{i}$,
respectively, and $D:=\mathsf{Diag}(s_{1},s_{2},\ldots,s_{N})$ be a diagonal
matrix containing the canonical correlations. We now state the main theorem:

######  Theorem 2.

Suppose $r_{t+1}$ and $x_{t}$ are $N$-dimensional jointly Gaussian random
variables. Then the objective function (2) is maximized at

$\displaystyle A=Q_{x}\mathsf{Diag}(L_{1},\ldots,L_{N})Q_{r}^{\prime},$ (46)

where diagonal matrix consist of the following elements

$\displaystyle L_{i}=\frac{1}{\gamma}\frac{s_{i}}{1+s_{i}^{2}},$ (47)

for $i=1,\ldots,N$.

Without loss of generality, let us reparameterize the matrix of coefficients
in terms of the canonical directions as $A=Q_{x}LQ_{r}^{\prime}$ where $L$ is
a variable matrix that we now have to optimize on. The expected value of the
portfolio returns at time $t$ can be written as

$\displaystyle\mathbb{E}[x^{\prime}_{t}Ar_{t+1}]$
$\displaystyle=\mathsf{Tr}(A\Sigma_{rx})$ (48)
$\displaystyle=\mathsf{Tr}(LD),$ (49)

and the variance as

$\displaystyle\mathsf{Var}[x_{t}^{\prime}Ar_{t+1}]$
$\displaystyle=\mathsf{Tr}(\Sigma_{x}A\Sigma_{r}A^{\prime})+\mathsf{Tr}(\Sigma_{rx}A\Sigma_{rx}A)$
(50) $\displaystyle=\mathsf{Tr}(LL^{\prime})+\mathsf{Tr}(DLDL).$ (51)

Putting all together, the objective function (4) can be written as:

$\displaystyle\max_{L}\mathsf{Tr}(LD)-\frac{\gamma}{2}(\mathsf{Tr}(LL^{\prime})-\mathsf{Tr}(DLDL)).$
(52)

Using the rules for matrix derivatives from Lütkepohl, (1997), we take first-
order conditions with respect to the matrix $L$ to get

$\displaystyle D-\gamma(L+DL^{\prime}D)=0.$ (53)

Since the matrix $D$ is diagonal and thus symmetric, we necessarily have the
following relationship

$\displaystyle\frac{1}{\gamma}D=L+DL^{\prime}D=L^{\prime}+DLD.$ (54)

Let $L^{a}:=L-L^{\prime}$, which is an anti-symmetric matrix (that is,
$L^{a^{\prime}}=-L^{a}$) and has diagonal elements zero. Then we have

$\displaystyle L^{a}=DL^{a}D.$ (55)

If we restrict our attention to the off-diagonal elements of (55), we see that
$L^{a}_{ij}=s_{i}s_{j}L^{a}_{ij}$ for all $i\neq j$, so
$(1-s_{i}s_{j})L^{a}_{ij}=0$. Therefore, either $s_{i}s_{j}=1$ for all $i\neq
j$ or $L^{a}_{ij}=0$. Consequently, $L^{a}$ is identically zero and hence, $L$
must be a symmetric matrix.

The benefit of variable matrix $L$ being symmetric is that the condition (53)
can now be written as

$\displaystyle L+DLD=\frac{1}{\gamma}D.$ (56)

If we focus on off-diagonal elements, we see that $(1+s_{i}s_{j})L_{ij}=0$ for
all $i\neq j$. But we also know that the canonical correlations satisfy the
following ordering $0\leq s_{1}\leq\ldots\leq s_{N}\leq 1$, and so it is
impossible that $s_{i}s_{j}=-1$. Thus, $L_{ij}=0$ for all $i\neq j$, and so
$L$ must be a diagonal matrix. Thus, optimizing over the elements
$L_{ii}\equiv L_{i}$ for $i=1,\ldots,N$, boils down to maximizing the
following univariate problems

$\displaystyle\max_{L_{1},\ldots,L_{N}}\sum_{i=1}^{N}L_{i}s_{i}-\frac{\gamma}{2}(L_{i}^{2}+s_{i}^{2}L_{i}^{2}).$
(57)

Our problem can now be easily solved to give us

$\displaystyle L_{i}=\frac{1}{\gamma}\frac{s_{i}}{1+s_{i}^{2}}.$ (58)

Hence, the diagonal elements are a nonlinear function of the canonical
correlations. The optimal matrix $A$ in (46) is effectively the optimal
scaling of the canonical portfolios of the asset returns and signals. ∎

  * Proof of Proposition 2.

Assuming that the second term in Equation (4) is approximately zero, we can
rewrite the objective function (5) in unconstrained form as

$\displaystyle\max_{\tilde{A}}\mathsf{Tr}(\tilde{A}\Sigma_{\tilde{r}\tilde{x}}^{\prime})-\frac{\gamma}{2}\mathsf{Tr}(\tilde{A}\tilde{A}^{\prime}),$
(59)

where we have used a change-of-variables
$\tilde{A}:=\Sigma^{1/2}_{x}A\Sigma^{1/2}_{r}$. The problem (59) can be solved
to yield

$\displaystyle\tilde{A}=\frac{1}{\gamma}\Sigma^{-1/2}_{x}\Sigma_{rx}^{\prime}\Sigma^{-1/2}_{r}.$
(60)

Since the solution is expressed in a different basis, we rescale it back to
the original assets through the following operation

$\displaystyle
A=\Sigma^{-1/2}_{x}\tilde{A}\Sigma^{-1/2}_{r}=\frac{1}{\gamma}\Sigma^{-1}_{x}\Sigma_{rx}^{\prime}\Sigma^{-1}_{r},$
(61)

As a result of this approximation, we now have

$\displaystyle L_{i}\approx\frac{1}{\gamma}s_{i}.$ (62)

Figure 2 contrasts the optimal singular value adjustment in (47) to the linear
approximation (62) for $\gamma=1$. We can see that for small values of
$s_{i}$, the optimal adjustment is approximately linear while for large values
of $s_{i}$, the optimal adjustment deflates the values of the canonical
correlations. ∎

  * Proof of Proposition 3.

We start by writing down the expression for the Lagrangian

$\displaystyle\mathsf{Tr}(A\Sigma_{rx})-\frac{\gamma}{2}\mathsf{Tr}(\Sigma_{x}A\Sigma_{r}A^{\prime})+\lambda(1-\mathbf{1}^{\prime}A^{\prime}x_{t}),$
(63)

$\lambda>0$ is the Lagrange multiplier. Using the change-of-variables
$\tilde{A}:=\Sigma^{1/2}_{x}A\Sigma^{1/2}_{r}$, we have

$\displaystyle\mathsf{Tr}(\tilde{A}\Sigma_{\tilde{r}\tilde{x}})-\frac{\gamma}{2}\mathsf{Tr}(\tilde{A}\tilde{A}^{\prime})+\lambda(1-\mathbf{1}^{\prime}\Sigma_{r}^{-1/2}\tilde{A}^{\prime}\Sigma_{x}^{-1/2}x_{t}).$
(64)

Performing the first-order conditions

$\displaystyle[\text{FOC }\tilde{A}]:$
$\displaystyle\enskip\Sigma_{\tilde{r}\tilde{x}}^{\prime}-\gamma\tilde{A}-\lambda\Sigma_{x}^{-1/2}x_{t}\mathbf{1}^{\prime}\Sigma_{r}^{-1/2}=0,$
(65) $\displaystyle[\text{FOC }\lambda]:$ $\displaystyle\enskip
1-\mathbf{1}^{\prime}w_{t}=0.$ (66)

By reverting the change-of-variables we have made and using the fact that the
portfolio policies are linear the signals, the portfolio weights is then

$\displaystyle w_{t}^{\mathrm{FI}}$
$\displaystyle=\frac{1}{\gamma}\Sigma_{r}^{-1}\Sigma_{rx}\Sigma_{x}^{-1}x_{t}+\lambda(x_{t}^{\prime}\Sigma^{-1}_{x}x_{t})\mathbf{1}^{\prime}\Sigma_{r}^{-1}.$
(67)

Solving for $\lambda$ gives

$\displaystyle\lambda=\frac{1-\kappa}{(x_{t}^{\prime}\Sigma^{-1}_{x}x_{t})(\mathbf{1}^{\prime}\Sigma_{r}^{-1}\mathbf{1})},$
(68)

where
$\kappa:=\gamma^{-1}\mathbf{1}^{\prime}\Sigma_{r}^{-1}\Sigma_{rx}\Sigma_{x}^{-1}x_{t}$.
Inserting the $\lambda$ into (67) and rearranging the expression gives us the
optimal portfolio (18). ∎

  * Proof of Proposition 4.

If $r_{t+1}$ and $x_{t}$ are $N$-dimensional zero-mean Gaussian random
variables, then their projections on to the $i$th canonical portfolios are
$u_{i}^{\prime}\tilde{r}_{t+1}$ and $v_{i}^{\prime}\tilde{x}_{t}$ are also
Gaussian random variables with mean zero and unit variances, respectively and
a joint correlation $s_{i}$. For the expected return of the $i$th canonical
portfolio, we have

$\displaystyle\mathbb{E}[\pi_{ti}]=s_{i}v_{i}^{\prime}\Sigma_{\tilde{r}\tilde{x}}^{\prime}u_{i}=s_{i}^{2}.$
(69)

Furthermore, the variance of the $i$th canonical portfolio is

$\displaystyle\mathsf{Var}[\pi_{ti}]$
$\displaystyle=\mathbb{E}[s_{i}^{2}(v_{i}^{\prime}\tilde{x}_{t})^{2}(u_{i}^{\prime}\tilde{r}_{t+1})^{2}]-s_{i}^{4}$
(70) $\displaystyle=s_{i}^{2}(1+2s_{i}^{2})-s_{i}^{4}$ (71)
$\displaystyle=s_{i}^{2}(1+s_{i}^{2}).$ (72)

The second equality follows from the second-moment-based result of the product
of two correlated Gaussian variables from Haldane, (1942, Section 6).
Substituting the expression of the optimal policy $A$ into the definition of
the expected value and variance of the portfolio returns, we have:

$\displaystyle\mathbb{E}[x_{t}^{\prime}Ar_{t+1}]=\frac{1}{\gamma}\mathsf{Tr}(\Sigma_{x}^{-1}\Sigma_{rx}^{\prime}\Sigma_{r}^{-1}\Sigma_{rx})$
$\displaystyle=\frac{1}{\gamma}\sum^{N}_{i=1}s_{i}^{2},$ (73)
$\displaystyle\mathsf{Var}[x_{t}^{\prime}Ar_{t+1}]=\frac{1}{\gamma}\mathsf{Tr}(\Sigma_{x}^{-1}\Sigma_{rx}^{\prime}\Sigma_{r}^{-1}\Sigma_{rx})$
$\displaystyle=\frac{1}{\gamma}\sum^{N}_{i=1}s_{i}^{2},$ (74)

where we assumed that the squared managed portfolio returns are approximately
zero in the expression for the variance term. ∎

  * Proof of Corollary 1.

This follows immediately Proposition 4 and the definition of the Sharpe ratio.
∎

  * Proof of Proposition 5.

Substituting the second-moment matrix of managed portfolio returns (38) into
the expression of the expected portfolio returns from Proposition 1 and using
the cyclic property of the trace operator, we have

$\displaystyle\mathsf{Tr}(\Sigma_{x}^{-1}\Sigma_{rx}^{\prime}\Sigma_{r}^{-1}\Sigma_{rx})=\frac{1}{\gamma}(\mu_{r}^{\prime}\Sigma_{r}^{-1}\mu_{r})(\mu_{x}^{\prime}\Sigma_{x}^{-1}\mu_{x})+\frac{1}{\gamma}\sum^{N}_{i=1}s_{i}^{2},$
(75)

and the claim follows. ∎