\section{Introduction}
\label{sec:Intro}
In the last five years, the Private Equity (PE) market saw a sharp expansion of its total economic value due to an unprecedent increase of PE investments  \citep{bainPEreport}.
However, predicting the outcome of these investments is a very challenging task, since
the future status of a private company can be extremely uncertain and it can range from bankruptcy, which corresponds to a total loss of the initial investment, to a Initial Public Offering (IPO), the most rewarding event for a PE investor. An impressive example is the investment of Peter Thiel in 2004: \$500,000 which brought him a 693,3\% return by the time of Facebook IPO \citep{facebook_worth}.
In addition to this extreme volatility of the investment outcome, PE investors have to struggle also with a severe lack of reliable and publicly available information related to the company
organizational and financial health.
For these reasons, PE fund managers and investors still strive to build  effective strategies
for  identifying private companies that will provide the highest investment returns.
Currently, most investors rely on naive methods, such as portfolio diversification, to reduce the investment risk.

In this paper, we tackle the PE exit prediction issue by proposing a data-analytic procedure that estimates the probability of a company to undertake an IPO in a given interval of time. The proposed method is based on an interplay of two models: a neural network (NN) that estimates the overall probability of the event of interest, and a Survival Model (SM) that more finely models this probability as a function of the forward time period.
Neural networks have been already used in the financial context, see, e.g.,  \citep{fadlalla2001analysis}, in particular they have shown to provide quite strong performance in exctracting valuable information from quantitative data \citep{kim2003discovery,samkin2008adding}.
Instead, Survival Models are broadly employed in the medical field \citep{james2013introduction,swan1996abstinence}, in particular for studying the probability of response to medical treatments in humans and animals, the development of diseases, and the identification and evaluation of prognostic factors and risks related to specific diseases, see \citep{elandt1980survival}.
SMs have been exploited also in many fields far from the medical ones, such as criminology, industrial production, insurance, and economics, see  \citep{james2013introduction,pitacco2004survival,modarres2016reliability}. To the best of the authors' knowledge, there is however no previous application of SMs to the PE market analysis; in particular,
among the survival approaches, the Accelerated Failure Time (AFT) model \citep{wei1992accelerated}
proved to be the most compatible with the PE framework of interest in this paper.
In a previous work \cite{calafiore2019classifiers} the authors explored a {\em static} model for a private company's exit prediction.
A time-to-IPO analysis, rather than a static analysis, adds an important tool to the PE firm toolbox, in as far as it allows to time appropriately potential intervention. Simply knowing if a potential target will experience an IPO or not, allows a PE firm to decide to invest or not. However, knowing that a company will have an IPO with a given probability in, say, one to two year time, gives to a PE firm a very real opportunity and timing to jump in and invest.

This paper discusses a time-to-IPO probabilistic analysis for PE firms. The experimental approach that we employed  is summarized as follows:
\begin{enumerate}
\item We considered a Thomson Reuters database containing relevant information  for a large number of companies, see Section~\ref{sec:data} for further details on the database and the features used for constructing the model.
\item Companies are labelled with 4 possible statuses: Bankrupt, Acquisition, IPO, Private. These classes synthesise the four main relevant  outcomes of a company in the PE market.
\item A neural network is trained to estimate the probability for each company of the union event of going bankrupt or being acquired.
\item We process the dataset by excluding bankruptcies and acquired companies. The dataset filtered in this way contains only companies that went IPO during the observation period or that are still private by the end of it, and constitutes our conditional dataset.
\item An Accelerated Failure Time (AFT) model is trained on the conditional dataset in order to model, for each company, the probability of occurrence in time of an IPO event. Since we fitted the AFT models on a dataset without bankruptcies and acquisitions, the probability obtained is conditioned on the knowledge that the company has not gone bankrupt or been acquired.
\item The total probability of going IPO in time is eventually computed by de-conditioning the AFT results.

\end{enumerate}

The contribution of this paper is the proposition, description and practical application of a machine learning procedure to the PE field.
We have tested this procedure on a set of real data, and we report
statistically significant results that may  have a high potential impact in this application field.

\section{Input Data}\label{sec:data}
The dataset used in this paper has been extracted from the financial platform Thomson Reuters Eikon \textsuperscript{\textregistered} \citep{eikonPage}.
This dataset contains the investment data of all the companies founded between 1998 and 2018 in Europe and in the United States. The records with missing data have been excluded.
After this procedure, we obtained a dataset of 47907 companies belonging to nine business sectors, as shown in Table~\ref{tab:exit_table}.
The available data for each company are: the names of the investors of the first 3 rounds of investment, the foundation date and the IPO date (if the company went public), the number of firms in each round, and the VIX index \citep{vixPage} value on the round dates, which is a public market sentiment indicator that measures the volatility (i.e., the uncertainty) of the public financial market.
Moreover, each company record is labeled with its industrial sector and the actual company status, a categorical value that has been aggregated in 4 labels: IPO, Bankrupt, Acquisition and Private.
The qualitative information of the investor names has been transformed into  quantitative data by means of an investor ranking approach: each investing firm is assigned to a ranking according to the number of rounds in which it has participated. The motivation behind investor ranking is simple: the largest the number of rounds an investor has performed, the strongest is its position in the PE market.
This approach has been employed also in other machine learning algorithms within the PE framework, see, e.g., \citep{bhat2011predicting,calafiore2019classifiers}.

\subsection{Descriptive statistics}
Figure~\ref{fig:Found_IPO_distr} shows the time distribution of IPOs and company foundations. It is worth to note how in the recent years the number of company foundations is substantially decreased with respect to the past. On the other hand, the number of IPOs is increased. This behavior shows the maturity of the company population considered in our analysis.
\begin{figure}
    \centering
    \includegraphics[width = 1\linewidth]{Figures/found_ipo_date_distr.png}
    \caption{Distribution of foundation (red line) and IPO dates (blue line) of the subjects (companies) under study.}
    \label{fig:Found_IPO_distr}
\end{figure}

Table~\ref{tab:exit_table} represents the final population distribution across status labels and business sectors. It is clear that our dataset suffers from  class unbalancing, since bankruptcies represents less than $2\%$ of the total data set. This issue has been handled with geometrical resampling techniques, as  described in Section~\ref{sec:NN}.

\begin{table}[]
\centering
\scalebox{0.8}{
\setlength\tabcolsep{3pt}
    \begin{tabular}{l|r|r|r|r|r|r}
\hline
 Sector & Bank. & IPO & LBO & M\&A & Private & \textbf{Total} \\
  \hline
  1 Communications & 165 & 1084 & 3 & 46 & 856 & 2154 \\
  2 Computer & 672 & 4926 & 4 & 140 & 11383 & 17125 \\
  3 Electronics & 151 & 940 & 1 & 27 & 1734 & 2853 \\
  4 Biotech & 84 & 1465 & 3 & 21 & 2704 & 4277 \\
  5 Medical & 50 & 1279 & 0 & 13 & 2599 & 3941 \\
  6 Energy & 5 & 492 & 1 & 9 & 852 & 1359 \\
  7 Consumer & 82 & 1851 & 5 & 33 & 3444 & 5415 \\
  8 Industrial & 108 & 2285 & 6 & 34 & 3755 & 6188 \\
  9 Other & 96 & 1605 & 4 & 24 & 2866 & 4595 \\
  \hline
  \textbf{Total} & 1413 & 15927 & 27 & 347 & 30193 & \textbf{47907} \\
  \hline
    \end{tabular}}
    \caption{Number of companies in each sector for each status within our dataset.}
     \label{tab:exit_table}
\end{table}

\section{Survival Analysis}
\label{sec:SurvCurv}
In this section, we present the basic concepts of Survival Analysis.
Survival Models attempt to estimate the time to a certain event in a probabilistic framework.
In our model, we are interested in estimating for each subject $i$ the ``survival time'' $T_i$, which is defined here as the period of time from the beginning of the subject observation (i.e., its foundation date) to the date $t_{E_i}$ of the event of interest $E_i$ (i.e., the company's IPO date).
In order to estimate the survival time $T_i$, SM can benefit from the knowledge of additional features regarding the characteristics of each subject. These subject's features will deeply depend on the application domain of the study. In this paper, we consider as additional features four measures of each investment round: the average, the maximum and the minimum ranking of investors who participated in that round, as well as the number of such investors. Moreover, exclusively for the first two rounds the VIX index value is also retained. In total, we consider 14 features for each company.

\subsection{Survival Models Assumptions}
\label{subsec:SM_assumptions}
SA is typically performed either  by means of  non-parametric Cox models \citep{cox1972regression}, or by  the parametric AFT models \citep{wei1992accelerated}.
The first category has the advantage that it does not need an underlying parametric probability distribution. However, Cox models rely on a rather strong  Proportional Hazard Assumption (PHA), which states that the ratio of the hazard functions of two subjects drawn from the population is a constant not depending on time \citep{cox1972regression}.
For each subject $i$, the hazard function $h_i(t)$ is the probability that $E_i$ will happen in the next, very short, time interval $(t,t+\delta t)$ \citep{lee2003statistical_haz}:
    \[
        h_i(t) = \lim_{\Delta t \to 0} \frac{ \mathbb{P}(t_{E_i}\  \in (t,t+\Delta t)|t_{E_i} \geq t) } {\Delta t}.
    \]
The PHA implies that
\[
   h_i(t) = h_0(t)g(\Vec{x_i}),
\]
where $\Vec{x_i}$ is the feature vector of subject $i$, $h_0(t)$ is a baseline hazard function that does not depend on the subject, and $g(\cdot)$ is a deterministic function of the subject features. The hazard function ratio of two companies is then
\[
   \frac{h_i(t)}{h_j(t)} = \frac{h_0(t)g(\Vec{x_i})}{h_0(t)g(\Vec{x_j})} = \frac{g(\Vec{x_i})}{g(\Vec{x_j})},

\]
where $\Vec{x_i}$ and $\Vec{x_j}$ are the feature vectors of  subjects $i,j$.
In the specific application considered in this paper this assumption is quite strong, requiring that the ratio of the probability of going IPO of two companies would depend just on the investors who invested in them, and not on time.

In order to verify if the considered dataset satisfies the PHA, we computed the Schoenfeld Residuals (SR) \citep{schoenfeld1982partial} for each covariate  and performed the statistical test (SR test) proposed by \cite{grambsch1994proportional}.
This test is based on the Schoenfeld residuals of a fitted Cox model weighted by its estimated covariance matrix. The formulation of the Schoenfeld residuals is not trivial, as it relies on the first derivative of the Cox regression log-likelihood function, and we refer the reader to \cite{schoenfeld1982partial} for additional details. Using the Schoenfeld Residuals, the PHA becomes equivalent to the hypothesis that the correlation between time and these residuals is null for all the covariates. Therefore, rejecting this hypothesis, even if only for one variable, equals to rejecting the possibility of using Cox models themselves.

The SR test has been performed on a
 business sector by sector basis.
In Table~\ref{tab:srtest_results} are reported the $p$-values of the SR test computed on all the covariates (Global $p$-values).

 For each sector, the Global $p$-value is well below $5\%$, therefore the PHA hypothesis can be globally rejected with a $95\%$ significance.
\begin{table*}
\centering
\scalebox{0.95}{
\setlength\tabcolsep{4.5pt}
    \begin{tabular}{c|c|c|c|c|c|c|c|c|c}
    \hline
 Sector & 1 & 2 & 3 & 4 & 5 &6 & 7 & 8 & 9 \\
  \hline
  Global $p$-value &1.21e-12 & 7.36e-72 & 3.88e-14 & 4.37e-08 & 1.20e-06& 4.76e-02 & 1.015e-11 & 5.47e-16 & 2.19e-07 \\
  \hline
    \end{tabular}}
    \caption{Results of the Schoenfeld residuals test on the Cox model as a whole for each sector.}
    \label{tab:srtest_results}
\end{table*}

\subsection{AFT model}
Since Cox model assumptions appeared to be invalidated by our data, we considered AFT models, which no not rely on such assumptions.
Within the AFT framework, each $T_i$ takes the form of a parametric random variable whose distribution is assumed to be known (e.g., Exponential, Weibull, lognormal etc.), while its parameters depend on the feature vector $\Vec{x_i}$, which is assumed to be directly related to the survival time $T_i$.

The AFT model uses a log-linear relation to describe the dependency between $T_i$ and the feature vector of the subject, see \citep{lee2003statistical_aft_theory}:
\begin{equation}
    \log T_i = a_0 + \sum^m_{j=1}a_{ij} x_{j} + \sigma\epsilon = \eta_i + \sigma\epsilon,
    \label{eq:aft_regr}
\end{equation}
where $x_{j}$, $a_{ij}$, and $a_0$ are respectively the $m$ covariates, the $m$ regression coefficients and the intercept parameter. All of these quantities can be summarized with $\eta_i$, which represents the deterministic term. On the other hand, $\epsilon$ is a stochastic error term, equipped with a density function $g(\epsilon)$, scaled by the unknown parameter $\sigma$. The error term $\epsilon$ is the actual source of randomness of the model and its probability distribution determines the resulting $T_i$ distribution.
It is important to underline that the regression coefficients and the parameter $\sigma$ are fundamental to compute the distribution parameters of each random variable $T_i$.
As an example, if $\epsilon$ has a standard normal distribution, then each $T_i$ will have a lognormal distribution having mean $\mu_i = a_0 + \sum^m_{j=1}a_{ij} x_{ij} = \eta_i$ and standard deviation $\sigma_i = \sigma$.
We compute the regression coefficients and the parameter $\sigma$ employing a Maximum Likelihood Estimation (MLE).
To test the statistical significance of the covariates, we used the Wald's statistic \citep{huber19721972}, whose results are reported in Section~\ref{sec:results}.

\section{Model fitting} \label{sec:ExpDesign}
In this section, we describe the procedure employed in order to fit the proposed survival model to the PE data.

\subsection{Data preprocessing}
When we apply SMs to PE data, the first problem encountered is the presence in the dataset of bankrupt and acquired companies. Indeed, plain SMs are able to model only one event and do not consider the case when the possible outcomes are more than one. More critically,
the timing information about bankrupt or acquisition events is rarely available in the database.
To deal with such multiple outcomes, we preprocess the dataset excluding the companies that went bankrupt or have been acquired. Then, we build the AFT model on the remaining observations (IPO and Private). This means that the probability that we estimate is conditioned by the fact that we know that the company has not experienced bankrupt or acquisition during the observation period. In addition, we divide the dataset according to the companies' industrial sector, and we train a different SM for each sector.

\subsection{Censoring}
A typical  problem that we face when we apply SMs to PE data is the presence of companies whose status remains  private within the data time span.
These companies enter the observation period at certain times (the foundation date) and by the end of the study the IPO event has not occurred yet.
We computed their survival times as the difference between the end of the observation period and their foundation date, and we marked with a binary variable these companies as censored observations.
The capability of dealing with censored data is actually a key property of Survival Models.
The censoring label encodes
the fact that for a censored company we do not know the exact survival time, but nevertheless we know that this survival time is {\em larger} than the observation time, and this information is exploited in the computation of the survival probability, see, e.g., \citep{lee2003statistical_censoring}.

\subsection{AFT procedure} \label{subsec:aft_procedure}
We use the AFT model to estimate the conditional survival probability $\mathbb{P}(T_i>t|V_i)$, where $V_i$ represents the event that the company $i$ does not experience bankrupt or acquisition during the observation period.
It is worth mentioning that in the analysis of the PE market the quantity of interest is the probability for a company $i$ of going IPO in a given time interval (i.e., $\mathbb{P}(T_i\leq t)$). However, since SMs are mainly used in the biomedical field, they are designed to model the survival probability $\mathbb{P}(T_i>t)$. Clearly,
it holds that
$\mathbb{P}(T_i\leq t) = 1 - \mathbb{P}(T_i>t)$.

In order to capture the different sector characteristics, we estimated a different SM for each of the nine industrial sectors. To train the nine AFT models, for each sector we have to made two choices: which distribution to choose for the stochastic error term $\epsilon$ and which covariate to use. We decided to test 4 different distributions in order to explore 3 different families: Exponential (Exponential, Weibull), Normal (LogNormal) and Snedecor's F (Generalized F) \citep{cox2008generalized}.

Our model selection procedure relies on a numerical (Likelihood, p-value)
and graphical (fitting plots) evaluation, which has been performed on a training dataset. First, for each distribution $k$ and for each sector $s$ a separate model $M_{k,s}$ has been fitted and the statistical significance of the covariates has been verified with the Wald test.
Then, each $M_{k,s}$ has been fitted again on the $90\%$ significant variables found previously for that model itself and its Maximum Log-Likelihood Estimate (MaxLLE) recorded.
Finally, the unconditioned AFT survival fitting curves have been plotted against the validation set curves (Figure~\ref{fig:aft_val_uncond}), that were obtained through a Kaplan-Meier \citep{kaplan1958nonparametric} empirical estimate (KM curve). The latter computation relies on a simple and robust counting process that is widely recognized as a consistent reference to benchmark graphical accuracy of the AFT fitted models \citep{lee2003statistical}. The results of this model selection procedure are discussed in Section~ \ref{sec:results}.

After having estimated the conditional probability $\mathbb{P}(T_i>t|V_i)$ with the procedure described above, we compute the marginal probability $\mathbb{P}(T_i\leq t)$.
By the law of total probability, the following holds:
\begin{equation}
\begin{split}
\mathbb{P}(T_i \leq t) &=1-\mathbb{P}(T_i > t) \\
&= 1- \mathbb{P}(T_i > t |V_i)
\mathbb{P}(V_i),
\end{split}
\label{eq:marg_prob}
\end{equation}
where $\mathbb{P}(V_i)$ has been estimated, for each company, using the neural network described in the following section. The second equality holds because we have assumed that $\mathbb{P}(T_i > t |\bar{V}_i)$=0, where $\bar{V}_i$ is the event that the $i$-th company has been acquired or has gone  bankrupt during the observation period. This seems a fairly reasonable assumption, since if a company was acquired or went bankrupt, it cannot go IPO, and the case where a public company go bankrupt is so rare that it is not worth taking account of it.

\subsection{Neural estimation}
\label{sec:NN}
In this section we describe the simple neural network (NN) approach that
was used to estimate the probability $\mathbb{P}(\bar{V}_i)$ that the $i$-th company is acquired or goes bankrupt. The input data for the NN are the features described in Section~\ref{sec:data} except for the Foundation date and the IPO date. Neural Networks have proven to achieve strong predictive performance within the financial framework \citep{fadlalla2001analysis}, especially for interpreting patterns from qualitative data \citep{kim2003discovery,samkin2008adding}. Moreover, it has already been shown that the selected features can constitute a strong basis for Machine Learning algorithms within the Private Equity framework \citep{bhat2011predicting}.
\begin{figure*}[tb]
    \centering

    \includegraphics[width = 1\textwidth]{Figures/aft_val_grande.pdf}
    \caption{Unconditioned AFT fitting compared to KM estimate built on the validation set, sector by sector. Black lines: Kaplan-Meyer empirical survival estimate of the validation set. Colored lines: best AFT fitted distributions on fitting test, for each sector. Weibull (blue line), Lognormal (green line), and Generalized F (yellow line) distribution. Dashed lines: $95\%$ confidence intervals.}
    \label{fig:aft_val_uncond}
\end{figure*}

The proposed neural network is a 1-layer Multilayer Perceptron (MLP). We use as activation function the Scaled Exponential Unit (SeLU) \citep{klambauer2017self}, since it provides strong evidences in preventing gradient vanishing optimization issues and promotes information flow through the net. Each company has been labeled as Bankrupt-Acquisition (BA) or Non-Bankrupt-Acquisition and the global dataset has been split into a training and a test set with a $2/3$--$1/3$ ratio.
Since the training dataset is  unbalanced  (the BA companies are 1787 out of 47907), a geometrical resampling technique has been used, namely the SVMSMOTE described in \citep{tang2008svms}.
The numerical optimization procedure has been performed using the Adam algorithm \citep{kingma2014adam}, a well-know method with self-adapting moment estimation, and the network has been trained for 70 epochs.

The performance metrics for both classes are shown in Table~\ref{tab:nn_results}. The Precision measures the ratio of elements correctly identified over the total elements predicted in that class, the Recall measures the ratio of elements actually identified over the total number of elements present in that class in the validation set, while the accuracy represents the overall precision measure in both classes. The output probability of the net was classified as {BA} if it was over $50\%$. The results obtained are comparable to the ones found in literature for a Machine Learning classifier on a Private Equity dataset \citep{bhat2011predicting,calafiore2019classifiers}, even if these related works act on a more balanced dataset.

\begin{table}[tb]
        \centering
        \scalebox{0.89}{
        \setlength\tabcolsep{4.5pt}
        \begin{tabular}{l|l|l|l|l}
        \hline
            Precision+ & Recall+ &  Precision- & Recall- & Accuracy\\
            \hline
            0.10 & 0.99 & 0.67 & 0.81 & 0.81  \\
            \hline
        \end{tabular}}
        \caption{NN predictive performance: ``+'' refers to BA event, ``-'' to Not-BA.
        }
        \label{tab:nn_results}
\end{table}

\section{Results}\label{sec:results}
In this section, we present the experimental results of the neuro-survival model described in the previous sections.
In order to fairly evaluate the proposed models, we divide the original dataset into a training and a validation set according to a $9/10$--$1/10$ proportion.

As discussed in Section~\ref{subsec:aft_procedure}, for each model $M_{k,s}$ we have to choose the distribution of $\epsilon$ and select the significant covariates. The Wald's test for each $M_{k,s}$ resulted in a number of significant covariates reported in
Table~\ref{tab:num_sel_covs_aft_completa}.
Instead, the distribution choice is the result of an integrated approach. It has to maximize the MaxLLE, minimize the number of Wald's selected covariates, and well approximate the KM curves.
As for the MaxLLE absolute value for the models fitted with the selected covariates, the difference across distributions in each sector was extremely low compared to its order of magnitude. So, this metric alone could not be considered as a robust approach for selecting the best fitting distribution. However, the distributions with the highest MaxLLE are the ones with the least number of selected variables in the Wald's test, in all sectors except sector 9. Then, in order to definitively select the distribution, we have analysed the fitting plots choosing the distributions  whose fitting curves  better capture the KM curves pattern. Table~\ref{tab:num_sel_covs_aft} reports the distribution and the number of covariates chosen for each sector.
The three criteria agree on the best distribution for all sectors except sectors 4 and 9.
In sector 4, the LogNormal distribution was found to approximate the KM curve better than the Genaralized F, which was the distribution with the highest MaxLLE. This was evident on the right tail of the KM curve, where the Genaralized F lied outside the KM confidence intervals, for this reason the LogNormal distribution was preferred.
In sector 9, according to the Wald's test, Exponential or LogNormal distribution were the best choices, however the Weibull distribution showed a better fitting and a higher MaxLLE. This partial contradiction can be related to the complex empirical pattern of this sector, since it groups all the companies without a specific industrial activity classification.

We tested the models obtained with this procedure on the validation dataset and computed their marginal probability $\mathbb{P}(T_i \leq   t)$ using Eq. \eqref{eq:marg_prob}. Figure~\ref{fig:aft_val_uncond} shows the fitting curves compared to the KM estimate. We can see that the obtained curves are strongly consistent with the validation data. In each sector, the selected models (colored lines) fit the empirical estimates (black lines)  reliably: fitting curves are within the KM estimates confidence interval, except for a few very restricted areas.

\begin{table}[tb]
\begin{center}
\begin{tabular}{l|rrrr}
  \hline
Sector & Exp. & Wei. & LogNorm. & Gen. F \\
  \hline
1 & 8 & 10 & 8 & 6 \\
  2 & 10 & 9 & 8 & 7 \\
  3 & 7 & 7 & 7 & 7 \\
  4 & 8 & 8 & 6 & 6 \\
  5 & 3 & 3 & 3 & 3 \\
  6 & 0 & 0 & 1 & 0 \\
  7 & 2 & 2 & 2 & 0 \\
  8 & 3 & 3 & 3 & 0 \\
  9 & 2 & 3 & 2 & 3 \\
   \hline
\end{tabular}
 \caption{Number of Wald's test $90\%$ significant covariates for each sector and distribution}
 \label{tab:num_sel_covs_aft_completa}
\end{center}
\end{table}

\begin{table*}[t]
\begin{center}
 \scalebox{0.83}{
 \setlength\tabcolsep{3pt}
\begin{tabular}{l|r|r|r|r|r|r|r|r|r}
  \hline
Sector (selected distr.) & 1 (Gen.F) & 2 (Gen.F) & 3 (Gen.F) & 4 (LogNorm.) & 5 (Gen.F) &6 (LogNorm.) & 7 (Wei.) & 8  (Wei.) & 9  (Wei.) \\
  \hline
\# sign. covariates &  6 & 7 & 7 & 6 & 3 & 1 & 2 & 3 & 3\\
   \hline
\end{tabular}}
\caption{Number of covariates with a statistical significance over 90\% and highest MaxLLE in that sector.}
\label{tab:num_sel_covs_aft}
\end{center}
\end{table*}

\section{Conclusions}
\label{sec:Conclusions}
In this paper, we presented a predictive model that estimates the time-to-IPO probability of private companies. The proposed method is based on an interplay of a neural network and a AFT survival model.
The results show that such neuro-survival model in the PE framework is able to reliably represent the probability in which a private company will go public in a given time interval.

The proposed method can be a useful support tool for PE investment decision and {\em timing}, enabling investors to estimate the probability of a company to experience a desirable IPO event within a given time interval or,
by complement, the probability that a PE investment will result in an adversary event.


\section*{Introduction}
\label{sec:intro}

Type theories with an impredicative universe \ensuremath{\mathsf{Prop}}\xspace of propositions,
such as the Calculus of Constructions \citep{coquandHuet:CoC},
lose the normalization property in the presence of a \emph{proof-irrelevant}
propositional equality $\_{=}\_ : \funT A \ensuremath{\mathsf{Type}}\xspace {A \to A \to \ensuremath{\mathsf{Prop}}\xspace}$
with the standard elimination principle.
The loss of normalization is facilitated already by a coercion
function with a reduction rule
\[
\begin{array}{l}
  \ensuremath{\mathsf{cast}}\xspace : \funT {A\,B} \ensuremath{\mathsf{Prop}}\xspace {\Eq \ensuremath{\mathsf{Prop}}\xspace A B \to A \to B} \\
  \cast A A e x \rhd x
\end{array}
\]
that does not inspect the equality proof $e : \Eq \ensuremath{\mathsf{Prop}}\xspace A A$ but only
checks whether the endpoints are (definitionally) equal.

The failure of normalization refutes a conjecture by Werner
\citeyearpar[Conjecture 3.14]{werner:lmcs08}.
Consistency and canonicity is not at stake; thus, the situation is
comparable to type theory with equality reflection
\citep{MartinLoef84,martinLoef:constrMathAndCompPrg}, aka
\emph{Extensional Type Theory}.
At the moment, it is unclear whether the use of impredicativity is essential to
break normalization; predicative type theory might be able to host a
proof-irrelevant propositional equality \citep{abel:nbe09}
while retaining normalization.

\section*{Counterexample to Normalization}

We employ the usual impredicative definition of absurdity $\bot$ and
negation $\neg A$ and a
derived definition of truth $\top$:
\[
\begin{array}{lcl}
  \bot   & : & \ensuremath{\mathsf{Prop}}\xspace \\
  \bot   & = & \funT A \ensuremath{\mathsf{Prop}}\xspace A
\end{array}
\qquad\qquad
\begin{array}{lcl}
  \neg   & : & \ensuremath{\mathsf{Prop}}\xspace \to \ensuremath{\mathsf{Prop}}\xspace \\
  \neg A & = & A \to \bot
\end{array}
\qquad\qquad
\begin{array}{lcl}
  \top   & : & \ensuremath{\mathsf{Prop}}\xspace \\
  \top   & = & \neg \bot \\
\end{array}
\]
The presence of $\ensuremath{\mathsf{cast}}\xspace$ allows us to define self-application under
the assumption that all propositions all equal.  The self-application
term $\omega$ refutes this assumption.
\[
\begin{array}{lcl}
  \delta & : & \top \\
  \delta\;z & = & z\;\top\;z
\end{array}
\qquad
\begin{array}{lcl}
  \omega & : & \neg \funT {A\,B} \ensuremath{\mathsf{Prop}}\xspace {\Eq \ensuremath{\mathsf{Prop}}\xspace A B} \\
  \omega\;h\;A & = & \ensuremath{\mathsf{cast}}\xspace\; \top\, A\; {(h\;\top\,A)}\;
    \delta
\end{array}
\]
Impredicativity is exploited in $\delta$ when applying $z : \bot$ to type
$\top = \bot \to \bot$ so that it can be applied to $z$ again.  The
type of
$\delta$
is $\top$ which we cast to $A$ thanks to
the assumption~$h$ that all propositions are equal.

We build a non-normalizing term $\Omega$ by applying $\omega$ to
itself through $\delta$, reminiscent of the shortest diverging term in untyped
$\lambda$-calculus.
\[
\begin{array}{lcl}
  \Omega & : & \neg \funT {A\,B} \ensuremath{\mathsf{Prop}}\xspace {\Eq \ensuremath{\mathsf{Prop}}\xspace A B} \\
  \Omega\;h & = & \delta\;(\omega\;h)

\end{array}
\]
Thanks to the reduction rule of $\ensuremath{\mathsf{cast}}\xspace$, term $\Omega\;h$ reduces to
itself:
\[
\begin{array}{lclcl}
  \Omega\;h
  & = & \delta\;(\omega\;h)
  & \rhd & \omega\;h\;\top\;(\omega\;h) \\
  & = & \cast \top \top {(h\;\top\;\top)} \delta \;(\omega\;h)

  & \rhd & \delta\;(\omega\;h) \\
  & = & \Omega\;h \\
\end{array}
\]
Thus, normalization is lost in the presence of a hypothesis (free
variable) $h$.  As a consequence, normalization that proceeds under
$\lambda$-abstraction can diverge.  This means that equality of open
terms cannot be decided just by normalization.

The counterexample can be implemented in Werner's
type theory with
proof-irrelevance \citeyearpar{werner:lmcs08}, refuting the normalization
conjecture (3.14).  We implement $\ensuremath{\mathsf{cast}}\xspace$ as instance of Werner's
more general equality elimination rule:
\[
\begin{array}{lcl}
  \ensuremath{\mathsf{Eq\_rec}}\xspace & : &
   \funT A \ensuremath{\mathsf{Type}}\xspace \funT P {A \to \ensuremath{\mathsf{Type}}\xspace}
    \funT {a\,b} A {P\;a \to \Eq A a b \to P\;b}
   \\
  \Eqrec A P a b x e & \rhd & x \qquad \mbox{if } a = b
\\[1ex]
  \cast A B e x & = & \Eqrec \ensuremath{\mathsf{Prop}}\xspace {(\lam a A \ensuremath{\mathsf{Prop}}\xspace)} A B x e
\end{array}
\]
The term $\Omega$ also serves as counterexample to normalization in
the theorem prover Lean \citep{deMoura:cade15},
version 3.4.2 \citep{lean342}.
\begin{adjustwidth}{\mathindent}{0em}
\lstinputlisting{propIrrEq.lean}

\end{adjustwidth}
Infinite reduction can now be triggered with the command
\lstinline{#reduce Omega}, which diverges.

\section*{A Counterexample Using Propositional Extensionality}

The counterexample of the last section used the absurd assumption that
\emph{all} proposition are equal.  The following counterexample
utilizes just the axiom of propositional extensionality, $\mathsf{propext}$,
which is a default axiom of Lean.
In fact, the weaker statement $\mathsf{tautext}$, which states the equality
of \emph{true} propositions, is sufficient.
\[
  \begin{array}{lcl}
\mathsf{propext} & : & \funT{A\,B} \ensuremath{\mathsf{Prop}}\xspace {(A \leftrightarrow B) \to \Eq \ensuremath{\mathsf{Prop}}\xspace A B} \\
\mathsf{tautext} & : & \funT{A\,B} \ensuremath{\mathsf{Prop}}\xspace {A \to B \to \Eq \ensuremath{\mathsf{Prop}}\xspace A B} \\
  \end{array}
\]
The counterexample uses the standard impredicative definition of
truth,
\[
  \begin{array}{lcl}
\top & = & \funT A \ensuremath{\mathsf{Prop}}\xspace {A \to A}
  \end{array}
\]
and a cast from
$\top \to \top$ to $A$, which are both tautologies under the
assumption $a : A$.
\[
\begin{array}{lcl}
\mathsf{id},\; \delta   & : & \top \to \top \\
\mathsf{id}\;x & = & x \\
\delta\;z & = & z\;(\top \to \top)\;\mathsf{id}\;z
\end{array}
\quad
\begin{array}{lcl}
\omega,\; \Omega & : & \top \\
\omega\;A\;a & = &
  \ensuremath{\mathsf{cast}}\xspace\; (\top \to \top)\; A\;
    (\mathsf{tautext}\;(\top \to \top)\;A\;\mathsf{id}\;a)\;
    \delta \\
\Omega & = & \delta\;\omega
\end{array}
\]
These definitions can be directly replayed in Lean 3.4.2 with the
standard prelude, yielding a
non-normalizing term \lstinline|Omega|.
\begin{adjustwidth}{\mathindent}{0em}
\lstinputlisting{propExtTrue.lean}
\end{adjustwidth}
Note that term \lstinline|Omega| is closed with respect the standard axioms
of Lean, and does not even have a weak head normal form.

\section*{Related Work and Conclusions}

The $\ensuremath{\mathsf{cast}}\xspace$ operator is inspired by Girard's operator
$J : \funT {A\,B} \ensuremath{\mathsf{Prop}}\xspace {A \to B}$ with reduction rule
$J\;A\;A\;M \rhd M$ that destroys the normalization property of
System~F \citep{girard:ssls70,harperMitchell:ipl99}.
In contrast to $J$, our $\ensuremath{\mathsf{cast}}\xspace$ also requires a proof of
equality of $A$ and $B$, but this proof is not inspected and thus does not
block reduction if it is non-canonical.  Thus, the simple lie that all
propositions are equal is sufficient to trigger divergence.

Historically, the Automath system AUT-4 is maybe the first type-theoretic proof
assistant to feature proof-irrelevant propositions
\citep{deBruijn:automath4}.  The terminology used by de Bruijn is
\emph{fourth degree identification}, where proofs are expressions
considered to have degree 4, propositions and values degree 3, types
and $\ensuremath{\mathsf{Prop}}\xspace$ degree 2, and the universe $\ensuremath{\mathsf{Type}}\xspace$ of types degree 1.

Lean's type theory \citep{carneiro:thesis} features an impredicative
universe of proof-irrelevant propositions which hosts both
propositional equality and the accessibility predicate \citep[1.2]{aczel:inductive75}.  As both may
be eliminated into computational universes, decidability of
definitional equality is lost, as
demonstrated by Carneiro \citeyearpar{carneiro:thesis}
for the case of accessibility.  As a consequence,
typing is not decidable.

The type-theoretic proof assistants Agda and Coq have recently
\citep{gilbertCockxSozeauTabareau:popl19}
been equipped with a proof-irrelevant universe of propositions
(``strict \ensuremath{\mathsf{Prop}}\xspace'').  In this universe, propositional equality can be
defined, but cannot be eliminated into types that are not strict
propositions themselves.  Under this restriction,
Gilbert \citeyearpar[4.3]{gilbert:PhD}
formally proved normalization and decidability of type
checking for the predicative case.

Several open problems remain:
\begin{enumerate}
\item Does the theory with impredicative strict \ensuremath{\mathsf{Prop}}\xspace have
  normalization and decidability of type checking as well?
\item Does the addition of Werner's rule, while destroying proof
  normalization, preserve decidability of conversion and type
  checking?  (Since proofs
  are irrelevant for equality, they need not be normalized during type
  checking.)
\item Does Werner's rule preserve normalization in the predicative
  case?  (Our counterexamples make use of impredicativity.)
\end{enumerate}

\subsection*{Acknowledgments}
\noindent
The authors acknowledge support by the Swedish Research Council
(Vetenskapsrådet) under grants 2014-04864 \emph{Termination
  Certificates for Dependently-Typed Programs and Proofs via
  Refinement Types} and
2017-04064 \emph{Syntax and Semantics of Univalent Type Theory}.
Our research group is part of the EU Cost Action
CA15123 The European research network on types for programming and
verification (EUTypes).
We thank Mario Carneiro for contributing the original Lean implementation
of the first counterexample.

\bibliographystyle{alpha}


\section{Introduction}

Solar rotation is one of the essential properties of the Sun as a star, but it is still to be understood clearly from the existing large volumes of observational data \citep{1985SoPh..100..141S,2000SoPh..191...47B,2014AJ....148...12X}. Generally speaking, traditional measurements of solar differential rotation in the surface atmospheric layers could be carried out mainly by using three methods: magnetic tracers, spectroscopic observations, and flux modulations \citep{1974ARA&A..12...47G,1984ARA&A..22..131H,2001ApJ...548L..87V}. Besides, the understanding of the interior rotation of the Sun varying with latitude and radius is also very important, and is usually determined by using the techniques developed through helioseismology \citep{1985Natur.317..591B,2003ARA&A..41..599T,2009LRSP....6....1H}. However, each of these methods has its own difficulties and limitations. That is to say, the high-precision determination of solar rotation periodicity is not an easy work as it depends on both the types of observational data and the analysis approaches.

Previous studies showed that there is a general consensus on the magnitude and form of the interior rotation as inferred from helioseismology \citep{1998MNRAS.298..543A}, the surface differential rotation rate as derived from sunspots \citep{2013SoPh..287..197J}, Doppler velocity measurements \citep{1990ApJ...351..309S}, as well as the magnetic activity features \citep{1993SoPh..145....1K,2016AJ....151...76X}, but there is no such consensus regarding the rotation rule for the features in the solar corona, as pointed out by \cite{2013ApJ...763..137H}. The main reasons why the coronal rotation is relatively less determined are as follows: i) there are virtually no obvious tracers, ii) the corona is the optically thin region, and iii) it is difficult to measure the coronal magnetic field \citep{1998SoPh..181..351V,2018MNRAS.475.3117B}. At present, it is widely accepted that the rotation of the solar corona could provide additional information on the construction of solar dynamo models, because it reflects the rotation characteristics in the sub-photospheric layers of the Sun \citep{2013IAUS..294..399K,2018AJ....156..152X,2018SoPh..293..128B}.

Lots of researchers have analyzed the coronal line emission at \normalsize{5303} \AA~ (Fe \tiny{XIV} \normalsize{green} line) to investigate the rotation characteristics of the solar corona \citep{1989ApJ...336..454S,1994SoPh..152..161R,1999SSRv...87..211I,2003SoPh..213...23A,2006AdSpR..38..906B,2010NewA...15..135B,2017MNRAS.466.4535B}. To reveal the rotational behavior of the solar corona in the Fe \tiny{XIV} \normalsize{green} line, \cite{1989ApJ...336..454S} analyzed the observational data made with the Sacramento Peak 40 cm coronagraph between 1973 and 1985. They found that the synodic rotation period is 27.52$\pm$0.42 days for the latitudes from 30N to 30S, and the corona rotates more rigidly than the features in the lower atmosphere. Using the homogeneous data set of coronal green line intensities over the period 1964-1990, \cite{1994SoPh..152..161R} calculated the averaged synodic rotation period as a function of latitude and time, and found that the coronal rotation period for the whole range of latitudes and for the latitudinal band $\pm30^\circ$ are 28.18$\pm$0.12 days and 27.65$\pm$0.13 days, respectively. Based on the data observed by the SOHO/LASCO C1 coronagraph during the solar minimum from April 1996 to March 1997, \cite{1999SSRv...87..211I} studied the periodicity and recurrence of Fe \tiny{XIV} \normalsize{emission} structures with heliospheric latitude and distance above the Sun's surface, they found no significant deviation from a rigidly rotating Fe \tiny{XIV} \normalsize{corona} with latitude or with distance from the Sun even on the small scales. Using the database of the coronal green line brightness (CGLB) covering six solar cycles (1939-2001), \cite{2006AdSpR..38..906B} found that the period of coronal rotation increases from 27 days at the solar equator to a little more than 29 days within the latitudes $\pm40^\circ$. Moreover, they found that, for the higher latitudes, the rotation period changes very slowly and its value of about 29 days remains up to the polar regions. That is to say, the total rotation of the solar corona could be described as a superposition of two rotation modes: the fast mode near the equator is around 27 days, and the slow mode exceeds 30 days \citep{2010NewA...15..135B}. Using the same but longer database, \cite{2017MNRAS.466.4535B} investigated the time variation of the differential rotation parameters of the solar corona, and found that the equatorial rotation rate of the corona increases in the epochs of minimum between the even and odd cycles and arrives at its minimum values between the odd and even cycles. They suggested that the 22-year cycle of differential rotation pattern could be used to explain the effect by the Gnevyshev-Ohl rule.

Meanwhile, solar radio fluxes at different frequencies are also used to determine the temporal variation of solar coronal rotation \citep{1997EM&P...76..141V,1998SoPh..181..351V,2001ApJ...548L..87V,2009MNRAS.400L..34C,2011MNRAS.413L..29V,2011MNRAS.414.3158C,2012MNRAS.423.3584L,2017ApJ...841...42X,2017SoPh..292...55B,2018MNRAS.475.3117B}. By using the disk-integrated simultaneous measurements of solar radio flux at the frequency range of 275-2800 MHz,  \cite{2001ApJ...548L..87V} found that the coronal rotation period strongly depends on the heights in the corona, i.e. the sidereal coronal rotation period increases with increasing frequency. However, by using the same data sets, \cite{2017SoPh..292...55B} found that the solar corona rotates more slowly at higher altitudes, which contradicts the conclusion given by \cite{2001ApJ...548L..87V}. By studying the coronal rotation period as a function of latitude between $\pm60^\circ$ for the period from 1999 to 2001, \cite{2009MNRAS.400L..34C} found that the solar corona rotates less differentially than the photosphere and chromosphere. The similar but abundant results were confirmed by \cite{2010MNRAS.407.1108C} who studied the differential rotation of soft X-ray corona derived from the solar full disk images of soft X-ray telescope on board the $Yohkoh$ observatory. Using the 10.7 cm solar radio emission data for the time interval 1947-2009, \cite{2011MNRAS.414.3158C} found that sidereal coronal rotation period varies from 19.0 days to 29.5 days with an average of 24.3 days, and a 22-year component (which might be related to the 22-year Hale magnetic cycle or magnetic field reversal) exists in the long-term rotation variation. Based on their analysis results,  \cite{2012MNRAS.423.3584L} investigated the long-term variation of the coronal rotation of radio emission, and found a weak decreasing trend in the coronal rotation period. They found that there is no 11-year Schwabe cycle of statistical significance for the secular variation of rotation period length, and the solar coronal rotation does not exhibit any systematic pattern. By applying the cross-correlation analysis and wavelet transformation, \cite{2017ApJ...841...42X} studied the daily and monthly radio flux at 2800 MHz from 1947-2014. They found that the coronal rotation period varies with the solar cycle phase, and the rotation period is relatively longer around the minimum year of a solar cycle. Interestingly, they found the coronal rotation variation is related to the 11-year Schwabe cycle, which is disagreement with the result obtained by \cite{2011MNRAS.414.3158C} and \cite{2012MNRAS.423.3584L}. In addition, the variability of fractal dimension of solar radio flux at different frequencies is also studied by \cite{1997EM&P...76..141V} and \cite{2018MNRAS.475.3117B}, and the main conclusion is that the fractal dimension increases with increasing frequency, i.e. randomness increases towards the inner corona.

Clearly, the long-term variation of solar coronal rotation and its relationship with the eleven-year solar activity cycle are not well understood, more studies on these topics are thus required and needed. With the hope to add further information on the temporal variation and the underlying processes of solar coronal rotation, we focus on modified coronal index for nearly eight solar cycles, introduced and provided by the Slovak Central Observatory in Hurbanovo, through several time series analysis approaches. A brief introduction of the observational data used in this paper is shown in the next Section. In Section 3 the statistical analysis results for coronal rotation determination are presented. And finally the conclusions and discussions are given in the last Section.

\section{Observational Data}

\subsection{Modified Coronal Index}
\label{sec:data}

The coronal green line (Fe \tiny{XIV} \normalsize{5303} \AA) is the most typical emission line from the inner corona of the Sun, and the regular observations of this prominent line have been processed since the year of 1939 \citep{2005JGRA..110.8106R}. By searching a suitable measurement from a Sun-observing space-based probe, \cite{2010SoPh..263...43L} introduced the modified coronal index (MCI) to replace the original coronal index, which was first proposed by \cite{1975BAICz..26..367R}.

As pointed out by \cite{2010SoPh..263...43L}, MCI represents the physical measure of solar cycle signature in the outer atmosphere and its advantage is the longest time sequence. This indicator describes the average irradiance of green line emitted in one steradian towards the Earth, and the daily values of MCI from 1939 January 1 to 2019 May 31 are shown in Figure 1. Here, the time series of MCI is downloaded from the website of the Slovak Central Observatory in Hurbanovo\footnote{http://www.suh.sk/obs/vysl/MCI.htm}, and it is expressed in ($10^{-6}\rm W/m^2$) units.

\begin{figure}
   \includegraphics[width=\columnwidth]{Fig1.eps}
    \caption{The daily values of modified coronal index (unit: $10^{-6}\rm W/m^2$) for the time interval from 1939 January 1 to 2019 May 31.}
    \label{fig:figure1}
\end{figure}

\begin{figure}
    \includegraphics[width=\columnwidth]{Fig2.eps}
    \caption{The histogram distribution of the probability density for the MCI measurements (the cyan blocks). The probability distribution can be described by the GEV distribution (the red line).}
    \label{fig:figure2}
\end{figure}

\begin{figure}
    \includegraphics[width=\columnwidth]{Fig3.eps}
    \caption{The cumulative probability distribution of the MCI measurements (the blue line) and the fitting values (the red line).}
    \label{fig:example_figure}
\end{figure}

\begin{figure}
   \includegraphics[width=\columnwidth]{Fig4.eps}
    \caption{The daily values of sunspot number for the time interval from 1939 January 1 to 2019 May 31.}
    \label{fig:figure1}
\end{figure}

\begin{figure}
    \includegraphics[width=\columnwidth]{Fig5.eps}
    \caption{The histogram distribution of the probability density for the DSN values (the cyan blocks). The probability distribution can be described by the GEV distribution (the red line).}
    \label{fig:figure2}
\end{figure}

\begin{figure}
    \includegraphics[width=\columnwidth]{Fig6.eps}
    \caption{The cumulative probability distribution of the DSN values (the blue line) and the fitting values (the red line).}
    \label{fig:example_figure}
\end{figure}

\subsection{Probability Distribution}
\label{sec:data}

Since we are working with the truncated sets, the truncated distributions on the fit --- building the data set from each probability distribution function (PDF) and cumulative distribution function (CDF) are used. To look for the probability distribution of the MCI values, they are divided into twelve bands with each band spanning 0.83. The numbers of MCI values within each band are counted, and then they are divided by their total numbers to describe the probability density.

Figure 2 displays the histogram distribution of the probability density of the MCI values. Here, the ordinate is the ratios defining the percentage probability of the MCI values for each band.

In order to find the best fit distribution of the data set, the distributions of the exponential, logistic, normal, and generalized extreme value (GEV) are examined. We found that, among the four distribution functions, the probability distribution of the MCI values can be described by the GEV distribution that models block maxima of the variables \citep{2007A&A...472..293A}.

 The red line shown in Figure 2 is the GEV distribution, and the most probable value of the MCI measurements is found to be 1.93. The cumulative probabilities of the MCI measurements (the blue line) and the fitting values (the red line) are shown in Figure 3, and one can see that the two curves match quite well.

 \subsection{Statistic Test of the Distribution}

Ultimately, we need to test the relative performance of the GEV distribution that is used to fit the data. To quantify the relative performance, the Kolmogorov-Smirnov (K-S) statistic test, which corresponds to the biggest difference between the observed and the empirical CDFs, is used.

\begin{equation}
\mathrm{KS} = \max \mid \mathrm{CDF}_{\mathrm{obs}}(x) -\mathrm{CDF}_{\mathrm{emp}}(x) \mid
\end{equation}

where $x_{\mathrm{trunc}} \leqslant x \leqslant \infty$, and $x_{\mathrm{trunc}}$ is the limit value below which data is not used in the fit.

Unfortunately, the GEV distribution does not pass the K-S test (in which the null hypothesis assumes that the MCI data is drawn by the GEV distribution), because the value of the KS (here is 0.0397) is larger than the Table of the percentage points of Kolmogorov statistics (the value for the 29371 data points is 0.0071 at the 95\% confidence level) that was provided by \cite{miller1956table}.  We also examined the other three distributions (exponential, logistic, and normal distribution), however, all of them did not pass the K-S test. That is, the distribution of the MCI data is not the GEV distribution, although the GEV distribution is the best one among the four distributions.

 \subsection{Comparison with the Sunspot Number Data}

 \cite{2002SoPh..206..401M} analyzed the correlation of the coronal index and other solar indices for the time period 1965-1997. They found that the coronal index could be used as a representative index of solar activity in order to be correlated with different periodic solar-terrestrial phenomena useful for space weather studies. To confirm the features that cause the coronal modulation are rooted in surface magnetic structures of the Sun, here we directly compared the distribution of the MCI data with the sunspot indicator.

The time series of the daily sunspot number (DSN) used in this work is freely downloaded from the World Data Center (WDC) --- Sunspot Index and Long-term Solar Observations (SILSO), Royal Observatory of Belgium, Brussels\footnote{http://www.sidc.be/silso/datafiles}. Since 2015 July 1, the original version of the DSN data has been replaced by a new entirely revised data set (the so-called version 2.0). For details, please refer to \cite{2014SSRv..186...35C} and \cite{2016SoPh..291.2629C}. Here, the time period from 1939 January 1 to 2019 May 31, the common time interval to the MCI data, is extracted. Figure 4 displays the daily values of sunspot number during the past 80 years (1939-2019).

Similar to the MCI data, we also divided the DSN time series into twelve bands, and then calculated the PDF for each band. Four distributions (exponential, logistic, normal, GEV) are used to find the best fit distribution of the DSN data, and the results are shown in Figures 5 and 6. The results indicate that the probability distribution of the DSN can be described by the GEV distribution, but it does not pass the K-S test. That is to say, the distribution of the DSN data is not the GEV distribution, although the GEV distribution is the best one among the four distributions.

From the above analyses, we can arrive at a conclusion that both the MCI and the DSN data show the similar distribution behavior, but they are not the GEV distribution.

However, it should be pointed out that the extreme value theory can be used to study the solar extreme events. For example, \cite{2017ApJ...839...98A,2018ApJ...853...80A,2019SoPh..294...67A} employed this technique to study the temporal variability of the sunspot number and the solar radio flux, and found that there is an upper bound for the distribution of the peak values of the photospheric and coronal indicators. Therefore, the coronal magnetic activity exhibits the similar kind of behavior as the sunspot activity, implying that the features that cause the coronal modulation are rooted in surface magnetic structures.

\section{Analysis Results}

\subsection{Period Length of Coronal Rotation}

To determine the coronal rotation period of MCI varying with the time, the codes of wavelet transform analysis provided by \cite{1998BAMS...79...61T}\footnote{http://atoc.colorado.edu/research/wavelets/} is utilized. Before the method is applied, the average values of MCI are normalized first, which is a procedure that subtracts the average value and then divides by the variance of the data set. In our analysis, the statistical significance test is carried out by assuming that the red noise exists in the given time series, which is also described in the above paper. After the normalized procedure, the resulting local power spectra of MCI are displayed in Figure 7. As we focus on studying the temporal variation of solar coronal rotation (around 27 days), the horizontal dashed black lines, representing the periods of 26 days and 38 days, are used to show the relative higher power belts.

\begin{figure}
    \includegraphics[width=\columnwidth]{Fig7.eps}
    \caption{Local wavelet power spectrum of the daily MCI, and the horizontal dashed black lines are the periods of 26 and 38 days. The dashed red line is the cone of influence (COI) where the edge affects might distort the pattern.}
    \label{fig:example_figure}
\end{figure}

As shown in Figure 7, it can be seen that the relative higher power belts appear around 27 days at the maximum times of solar cycles, as the signal-to-noise is relative higher when the MCI itself is higher. For the minimum times of solar cycles, the solar corona is also rotating. However, the signal-to-noise is relative lower at the periods of lower activity, so the rotational behavior is difficult to seen from the local wavelet power spectra. According to Figure 7, the global wavelet power spectra of daily MCI with the periodic scale up to 128 days are calculated and shown in Figure 8. From this figure one can clearly see that the coronal rotation period with a value of 27.5 days is the only one peak at the periodic scales shorter than 64 days, which is statistically significant above the 95\% confidence level (dashed red line). It is noteworthy that the periodicity obtained here is the synodic coronal rotation period ($\rm P_{synodic}$), so the sidereal coronal rotation period ($\rm P_{sidereal}$) is determined to be 25.57 days. The approximate conversion relation between them could be described by the following expression \citep{1995SoPh..159..393R,1996SoPh..168..211W,2014SoPh..289.1471S}:

\begin{equation}
\rm P_{sidereal} = \frac {365.26 \times P_{synodic}} {365.26+P_{synodic}}
\end{equation}

The synodic coronal rotation period (27.5 days) obtained by us is good agreement with the results reported by previous researchers, whose studies focused on the Fe \tiny{XIV} \normalsize{green} line for the period 1973-1985 \cite[27.52$\pm$0.42 days]{1989ApJ...336..454S}, the CGLB over the period 1964-1990 \cite[27.65$\pm$0.13 days]{1994SoPh..152..161R}, and the Fe \tiny{XIV} \normalsize{emission} structures observed by SOHO/LASCO C1 coronagraph \cite[27.5$\pm$1 days]{1999SSRv...87..211I}. However, the obtained sidereal coronal rotation period in this work (25.57 days) is slightly greater than the result (24.3$\pm$0.2 days) given by \cite{2011MNRAS.414.3158C} who studied the daily radio emission at 2.8 GHz between 1947 and 2009.

\begin{figure}
    \includegraphics[width=\columnwidth]{Fig8.eps}
    \caption{Global wavelet power spectrum of the daily MCI with the periodic scale up to 128 days, and the dashed red line is the 95\% confidence level.}
    \label{fig:example_figure}
\end{figure}

\begin{figure}
    \includegraphics[width=\columnwidth]{Fig9.eps}
    \caption{The temporal variation of the period length of the coronal rotation (PLCR) during the considered time interval. The bold solid red line shown in this figure is the linear regression line, and the dashed red lines are the 95\% confidence bounds.}
    \label{fig:example_figure}
\end{figure}

\subsection{Temporal Variation of Coronal Rotation}

Since the local wavelet power spectra provide detailed contents in the time-frequency space, the local periodicity or frequency of daily MCI at a particular time point could be estimated. That is to say, in the local wavelet power spectra shown in Figure 7, the rotation period length at a certain time point has the highest power spectrum among the periodic scales of 24-30 days, where are the length range of coronal rotation period located \citep{1984ApJ...283..373H,2011MNRAS.414.3158C}. Based on this conception, the rotation period length of sunspot numbers \citep{2011Ap&SS.331..441L,2012RAA....12..187X}, sunspot areas \citep{2011ApJ...730...49L}, solar mean magnetic field \citep{2017AJ....153..171X}, and 10.7 cm solar radio flux \citep{2017ApJ...841...42X} are studied and estimated. In this sense, the similar procedure could be also applied to obtain the period length of the coronal rotation (PLCR) of daily MCI at each time point during the considered time interval, and the extracted results are shown in Figure 9. The time series of the PLCR is smoothed with a one-year sliding window, and the rotation period shown here is also the synodic coronal rotation period (the average value is 26.69$\pm$0.68 days).

As Figure 9 displays, the PLCR varies from cycle to cycle, and a general trend with decrease in length is clearly seen. To look into the secular trend of the PLCR time series, we perform a linear regression analysis to fit the data set, and the regression line (indicated by the bold solid red line) is also shown in this figure. The dashed red lines indicate the 95\% confidence bounds for the linear regression analysis. The relation between the two parameters can be described as $p(t)$ = -0.0110 ($\pm$0.0003) $*$ $t$ + 27.14 ($\pm$0.02), where $t$ is the time (in year, here the start time is set to 0) and $p(t)$ is the PLCR varying with the $t$. From the linear regression analysis, the PLCR is found to have a general trend, which exhibits a linear decrease from 27.14 ($\pm$0.02) days to 26.25 ($\pm$0.01) days (the period difference is 0.89 ($\pm$0.03) days with a slope of about -0.0110 ($\pm$0.0003) day/year).

To quantitatively judge the result of fit, the goodness-of-fit with 95\% confidence is used. The parameter of Adjusted R-square (Adj. R$^2$) value is calculated and represents the most useful measure of the success of the statistical model. In practice, if this value is close to zero, the fit result can be considered as a good fit. Another statistical parameter, the RMSE (root-mean-squared error), which is known as the fit standard error and the standard error of the regression, is also used. Similar to the Adj. R$^2$ value, the RMSE value close to zero indicates a good fit. Here, the Adj. R$^2$ value and the RMSE are 0.1416 and 0.6316, respectively. These two parameters are not so large that the fit result can be considered as a good fit. Therefore, at least from a global point of view, the solar corona rotates more and more speedily, i.e. accelerates its rotation rate, during the studied time interval (1939-2019).

It is clear from Figure 9 that the long-term decreasing trend of PLRC time series we obtained is similar to those obtained by \cite{2009A&A...497..835H}, \cite{2011MNRAS.414.3158C},  \cite{2012MNRAS.423.3584L}, and \cite{2017ApJ...841...42X}. Interestingly, all of these studies used the daily 10.7 cm solar radio flux to reveal the rotational behavior of the solar corona. Since both the 10.7 cm radio flux and the Fe \tiny{XIV} \normalsize{green} line emission reflect the global characteristic of the solar corona, the long-term temporal variation of PLCR obtained by these two indicators are similar to each other. A similar inference could be drawn at the photospheric magnetic activity estimated by \cite{2011ApJ...730...49L} who studied the rotational cycle length of daily sunspot areas from 1849 to 2010. They found that the rotation length of the photospheric rotation have a secular trend showing a linear decrease by about 0.47 days. However, the period difference of the PLCR obtained by us (0.89 $\pm$0.03 days) is almost twice as large as that of \cite{2011ApJ...730...49L}. The possible reason for the significant difference of the two works is that the studied time interval differs from each other. The time interval of the sunspot number used in \cite{2011ApJ...730...49L} is from 1849 January 1 to 2010 February 28 (about 162 years), whereas the time interval of the MCI data used in work is only 81 years. By using the similar fit method, a smaller period difference can be seen from longer time series. This reason can be confirmed by pervious works. For example, the time interval of the sunspot number used in \cite{2009A&A...497..835H} is from 1849-2004 (about 156 years), their work gave a value of 0.66 days. That is, the period differences obtained by different authors are not caused by the uncertainties of the rotation period. So from the present work and earlier studies, we infer that sunspot activity, coronal flux, and green corona exhibit the similar global rotation behavior.

\subsection{Periodicity in the PLCR Time Series}

There are many time-frequency analysis methods to search the periodic variation of solar time series. One of such statistical techniques is the auto-correlation analysis, which has successfully been employed for determining the periodicities of solar activity indicators \citep[see][]{1969SoPh...10..135H,1998SoPh..181..351V,2001ApJ...548L..87V,2010MNRAS.407.1108C}. By applying this technique to the PLCR time series, the auto-correlation coefficients are plotted against the phase lags (in days), and the analysis result is shown in Figure 10. We can easily find that, for the phase shifts smaller than 4500 days, there are four prominent peaks when the phase shifts are 3.25, 6.13, 9.53, and 11.13 years. The dashed red line shown in Figure 10 is the 95\% confidence level line, so the four periodicities obtained by us are statistical significant.

The periodicity of 3.25 years is the so-called quasi-biennial oscillation (QBO; \citealt{2014SSRv..186..359B}), and the similar periods of 2-3 years were also found by \cite{1995SoPh..158..173J}, \cite{1999SoPh..184...41J}, \cite{2009SoPh..257...61J}, and \cite{2017ApJ...841...42X}. These studied focused on the temporal variations of solar rotation derived from the sunspot group data, Doppler-velocity measurements, and daily 10.7 cm radio flux. Our analysis results further support the existence of QBO in the periodic variation of solar coronal rotation. Although the physical mechanism of solar QBO has not been fully understood, it is believed to be intrinsic to the internal dynamo process. Moreover, previous studies showed that the source of the QBO of the solar magnetic field should be situated at the base of the solar convection zone. Most recently, \cite{2010ApJ...724L..95Z}, \cite{2015NatCo...6E6491M}, and \cite{2018ApJ...856...32Z} concluded that the solar QBO of magnetic activities should be caused by the magnetic Rossby wave instabilities. Therefore, it is an evidence that the rotational behavior of the solar corona shows the similar temporal variation as the lower atmospheric layers, and both of them might be connected with dynamical process at the base of the solar convection zone.

\begin{figure}
    \includegraphics[width=\columnwidth]{Fig10.eps}
    \caption{The auto-correlation coefficient of the daily PLCR time series for the determination of the existing periodicities. The dashed red line is the 95\% confidence level line for the statistical significant test.}
    \label{fig:example_figure}
\end{figure}

\begin{figure}
    \includegraphics[width=\columnwidth]{Fig11.eps}
    \caption{The cross-correlation coefficient between the monthly PLCR and the monthly MCI. The dashed red lines are the 95\% confidence level lines for the statistical significant test.}
    \label{fig:example_figure}
\end{figure}

The periodicity of 6.13 years in the PLCR time series is close to the periodicity of 6.1 years in the differential rotation parameters of Greenwich sunspot groups  \citep{1995SoPh..158..173J}. In the statistical study of \cite{2016SoPh..291.3485J}, a periodicity around 5.4 years was also revealed in the equatorial rotation rate of sunspot groups. They thought that the periodicity of 5.4 years should  be the 3.5 times of the synodic period (1.597 year) of Venus and Earth \citep{2013PRP.....1..147W}. Most recently. by studying the temporal variation of coronal rotation in the 10.7 cm radio flux, \cite{2017ApJ...841...42X} found a periodicity of 6.6 years, which was interpreted as the third harmonic of the 22-year magnetic cycle. A similar periodicity of 6.8 years was found by \cite{2018ApJ...855...84X} who studied the equatorial rotation rate of solar magnetic fields from 1975 September to 2008 April. If the periodicity around 6 years is indeed the 22-year magnetic cycle, then our result is in good agreement with the results given by \cite{2000JApA...21..167J}, \cite{2011MNRAS.414.3158C}, and \cite{2017MNRAS.466.4535B}. That is, the rotational behavior of the solar corona might be related to the 22-year Hale cycle or the magnetic activity reversal. Further studies are needed to reveal the physical origin of the periodicity around 6 years.

The existence of the quasi-eleven-year periodicity in the temporal variation of solar rotation is an open question. According to the statistical analyses of solar radio flux at 2.8 GHz, \cite{2011MNRAS.414.3158C} found that the temporal variation of rotation period does not have any systematic periodicity around 11 years, and there is either not or weak correlation between coronal rotation and solar activity. Later, \cite{2012MNRAS.423.3584L} also studied the cycle-related variation of coronal rotation period, and they claimed that there is no 11-year periodicity for the secular variation of rotation period length. However, \cite{2017ApJ...841...42X} restudied the temporal variation of coronal rotation based on the continuous wavelet transform and auto-correlation analysis, they thought that the long-term variation of coronal rotation period should be related to the 11-year Schwabe cycle, because a periodicity of 10.3 years is detected from the smoothed coronal rotation period. By studying the rotation characteristics of historical solar observations, including the Greenwich sunspot groups, the spectroscopic velocity data, the large-scale magnetic fields, and the photospheric magnetic maps, the quasi-eleven-year periodicity of solar rotation variation has also been reported by \cite{1997SoPh..170..389J}, \cite{1999SoPh..184...41J}, \cite{2001SoPh..201....1O}, \cite{2006SoPh..237..365B}, and \cite{2018ApJ...855...84X}. In this work, we find that the quasi-11-year periodicity (here are 9.53 and 11.13 years) indeed exists in the PLCR time series. So from this, our analysis results strongly support the standpoint that the temporal variation of coronal rotation has a close connection with the 11-year Schwabe cycle.

\subsection{Phase Relationship between Coronal Rotation and Solar Activity}

Earlier studies \citep[see][]{2001SSRv...95..227R,2002SoPh..206..401M,2005AdSpR..35..410M,2007SoPh..241..269M,2014SSRv..186..105E} have shown that the intensity of the green corona is highly correlated with other solar activity indicators as all of them are intrinsically inter-linked through solar magnetism. So, the MCI is applied here as an indicator to describe the temporal variability of solar magnetic activity. To better understand the phase relationship between coronal rotation and solar activity, the cross-correlation analysis between the monthly PLCR and the monthly MCI is performed.  Here, we used the monthly data set to remove many small peaks existing in the daily data set. Figure 11 displays the resulting cross-correlation coefficients varying with the relative phase lags (in months), in which the abscissa is the phase shift of the PLCR versus the MCI with forward shifts given positive values.

As Figure 11 shows, three significant peaks, above the95\% confidence level, within the time shifts of $\pm$132 months ($\pm$ 11 years) are clearly seen. When the time shifts are -60 months (-5.0 years), 19 months (1.6 years), and 101 months (8.4 years), the values of the cross-correlation coefficient reaches a local maximum. Therefore, the phase relationship between the coronal rotation and the solar activity is very complex, without any systematic pattern.

It is well known that the wavelet transformation analysis is a statistical tool widely used to study the nonlinear and non-stationary signals in the time-frequency space. One of its extension, the cross-wavelet transform (CWT), can be applied to determine the relative phase difference between the two time series \citep{2004NPGeo..11..561G}. This technique has been successful used in a broad range of research aspects of solar physics fields, such as \cite{2015AJ....150..171X}, \cite{2017ApJ...850..120R}, \cite{2017ApJ...851..141X}, \cite{2018JGRA..123.6148G}, \cite{2018ApJ...857...28K}, and so on.

\begin{figure*}
    \includegraphics[width=2.0\columnwidth]{Fig12.eps}
    \caption{Panel a): the CWT spectra of the monthly PLCR with the monthly MCI, and the dashed white line is the cone of influence (COI) where the wavelet transform suffers from the edges effects. The relative phase relationship is shown as the arrows with in-phase pointing right, anti-phase pointing left, and the PLCR leading the MCI by 90$^\circ$ pointing down. Panel b): the relative phase difference (in
months) between the monthly PLCR and the monthly MCI varying with the time. Panel c): the relative phase differences (in months) between the monthly PLCR and the monthly MCI as a function of the period scales.}
    \label{fig:example_figure}
\end{figure*}

In this work, the codes provided by \cite{2004NPGeo..11..561G}, available in GitHub\footnote{https://github.com/grinsted/wavelet-coherence}, are employed to examine the phase differences varies with the time and the period scale between the monthly PLCR with the monthly MCI. The CWT spectra between them are shown in the panel a) of Figure 12, in which the dashed white line showing the cone of influence (COI) where the local spectrum suffers from the edges effects. In this figure, the relative phase difference is described as the arrows with in-phase pointing right, anti-phase pointing left, and the PLCR lagging behind the MCI by 90$^\circ$ pointing up. From this figure, one can easily see that in the whole frequency ranges (both the high-frequency and the low-frequency components), the power spectrum displays the strong phase mixing, that is, the distribution of the arrow directions (i.e. phase angles) is almost random, implying no regular fluctuation pattern between the two time series.

Based on the local phase relations displayed in the panel a) of Figure 12, we calculated separately the relative phase differences as a function of the time point and the period scale. The results are shown in the panel b) and the panel c) of Figure 12, respectively. Here, the phase difference at a certain time point is simply calculated by averaging all phase differences at that time point in the whole periodical scales. This concept is also used to calculate the relative phase differences varying with the periodic scales. As the two panels display,  the phase relationship between the coronal rotation period and the solar magnetic activity is not only time-dependent but also frequency-dependent.

As shown in Figure 12, when the signal is weak (mostly at high-frequency components), the distribution of the phase is random. If the signal is strong (around the 11-year cycle band), there is a systematic trend in the phase. In the panel a) of Figure 12, the ranges of the 11-year cycle component are indicated as the dashed bold red lines. It is clearly seen that most of the arrows (expect for the time interval 1960-1972) point down, indicating that the PLCR time series leads the MCI time series. According to our results, we should be careful in choosing the suitable periodic ranges and the time interval when we need to study their phase relationship in detail. A similar conclusion has been obtained by \cite{2007A&A...475L..33D} who studied the scale-resolved phase coherence of hemispheric sunspot areas during the time interval from 1874 to 2007. They found that the coherent phase variables of hemispheric sunspot activity are found to exist for a small frequency band with periodicities around the dominating 11-year cycle. So from our analysis, it is reasonably inferred that the small mismatch in the dominant 11-year cycle band brings out the phase to drift.

\section{Conclusions and Discussions}

In the present study, we investigated the long-term temporal variations of solar coronal rotation with the use of the MCI during the time interval between 1939 January 1 and 2019 May 31. Firstly, the continuous wavelet transform was applied to calculate the coronal rotation period from a global point of view. Based on the conception proposed by earlier studies, the time series of the PLCR was extracted. Then, the secular decreasing trend and the possible periodicities of the PLCR time series are obtained by the linear regression analysis and the auto-correlation analysis, respectively. And finally, the cross-correlation analysis and the cross-wavelet transform were used to examine the complex phase relationship between the monthly PLCR and the monthly MCI.

From a global point of view, the synodic coronal rotation period obtained by us is good agreement with the results reported by previous researchers, whose studies focused on the coronal green line emission, such as \cite{1989ApJ...336..454S}, \cite{1994SoPh..152..161R}, and \cite{1999SSRv...87..211I}. However, the obtained value is slightly greater than the result given by \cite{2011MNRAS.414.3158C} and \cite{2012MNRAS.423.3584L} who studied the periodicities of the coronal rotation derived from the daily radio flux at 2800 MHz. The slightly difference of the rotation period length between the Fe \tiny{XIV} \normalsize{green} line emission and the solar radio flux at 2800 MHz can be interpreted as the following reason. The 10.7 cm solar radio flux is directly related to the total amount of the magnetic flux, which should be originated from the lower (inner) corona at the 60 000 km above the surface atmosphere \citep{2001ApJ...548L..87V,2017SoPh..292...55B}. However, the Fe \tiny{XIV} \normalsize{green} line emission arises in coronal material with a characteristic temperature of about 1.8 $\times 10^{6}$ K \citep{1969MNRAS.142..501J}, which is the typical temperature of most of the corona outside of the coronal holes. Besides, this temperature is also found over active regions, where the Fe \tiny{XIV} \normalsize{emission} reaches its highest intensity. That is, the Fe \tiny{XIV} \normalsize{green} line emission is considered to be related to both the active regions and the quiet corona \citep{1989ApJ...336..454S,2018AstL...44..727B}.

Previous studies \citep[see][]{2009A&A...497..835H,2011Ap&SS.331..441L,2012RAA....12..187X,2017ApJ...841...42X} on the secular trend of the solar rotation period have often revealed a general decreasing trend in the solar rotation period. It should be pointed out that these studies focused on the sunspot activity and the corona radio flux. Our analysis results, based on the database of the Fe \tiny{XIV} \normalsize{green} line emission, somewhat confirm and enhance the earlier findings. So from this work and earlier studies, we infer that the sunspot activity, the coronal radio flux, and the green corona exhibit the similar decrease trend in the rotation period over the past several decades, i.e., the Sun accelerates its global rotation rate in the long run. However, an opposite long-term trend (namely, a secular deceleration) for the solar rotation rate was also obtained by some authors. For example, \cite{2005SoPh..232...25J,2005ApJ...626..579J}, \cite{2004HvaOB..28...55B}, and \cite{2006SoPh..237..365B} found that the solar rotation rate exhibits a secular deceleration trend. The reason for different opinions might be possibly related to the different data analysis methods those are used. The time-frequency signal analysis is applied to estimate the solar rotation period by \cite{2009A&A...497..835H}, \cite{2011ApJ...730...49L}, \cite{2017ApJ...841...42X}, and this work, which is mainly caused by the solar large-scale active regions or active longitudes. However, the solar rotation rate studied by \cite{2005SoPh..232...25J,2005ApJ...626..579J}, and \cite{2004HvaOB..28...55B} is estimated from the position variation of the surface tracers (sunspot or sunspot group). The key problem of interpreting the analysis results derived from different tracers is whether the observed changes represent the global variation of the rotation or might be caused by the specific property of the used tracer, as previously pointed out by \cite{2006SoPh..237..365B}. Further studies are needed to reveal this problem in the future.

Four significant periodicities with the values of 3.25, 6.13, 9.53, and 11.13 years exist in the temporal variation of the coronal rotation period. The first one might be related to the solar QBO, which should be caused by the magnetic Rossby wave instabilities in the solar tachocline \citep{2010ApJ...724L..95Z}. The possible physical origin of the second periodicity (6.13 years) has not fully understood, but it can be interpreted as the third harmonic of the 22-year magnetic cycle. Our analysis result supports the standpoint that the rotational behavior of the solar corona might be related to the 22-year Hale cycle or the magnetic activity reversal, which has been found by \cite{2005ApJ...626..579J}, \cite{2011MNRAS.414.3158C}, and \cite{2017MNRAS.466.4535B}. During the past several years, the existence of the quasi-eleven-year periodicity in the temporal variation of solar rotation is an open question, and many authors have given differing opinions. Based on our analyses, the last two periodicities (9.53 years and 11.13 years) can be considered as the quasi-11-year Schwabe cycle. That is, the temporal variation of the solar coronal rotation has a close connection with the eleven-year Schwabe cycle.

The relative phase relationship between the PLCR and the MCI, obtained by the cross-correlation analysis, does not display any regular pattern. That is, their phase relationship is very complex, which is also confirmed and enhanced by the CWT analysis. For both the high-frequency and the low-frequency components, the power spectrum shows the strong phase mixing and the distribution of the arrow directions is almost random. By calculating the relative phase differences as a function of the time point and the period scale, we found that the phase relationship between the coronal rotation period and the solar magnetic activity is not only time-dependent but also frequency-dependent. It should be noticed that the suitable periodic component and the time interval should be carefully chose and considered. For a small range around the 11-year cycle band, there is a systematic trend in the phase, and the small mismatch in this band brings out the phase to drift. To better understand the relationship between the coronal rotation and the solar magnetic activity, it will be very interesting to further investigate their complicated phase relationship in the future.

\section*{Acknowledgements}

The database  of daily and monthly modified coronal index used in this study was introduced by B. Luk\' a\u c and M. Rybansk\' y, and can be downloaded from the website of the Slovak Central Observatory in Hurbanovo (http://www.suh.sk/obs/vysl/MCI.htm). The sunspot number data was obtained from the WDC-SILSO, Royal Observatory of Belgium, Brussels (http://www.sidc.be/silso/). This work is supported by the National Key Research and Development Program of China (2018YFA0404603), the Joint Research Fund in Astronomy (Nos. U1631129, U1831204, U1931141) under cooperative agreement between the National Natural Science Foundation of China (NSFC) and Chinese Academy of Sciences (CAS), the National Natural Science Foundation of China (Nos. 11873089, 11903009), the Youth Innovation Promotion Association CAS, the Yunnan Key Research and Development Program (2018IA054), the open research program of CAS Key Laboratory of Solar Activity (Nos. KLSA201807), the Key Laboratory of Geospace Environment, CAS, University of Science \& Technology of China, the major scientific research project of Guangdong regular institutions of higher learning (2017KZDXM062), and the youth project of science and technology research program of Chongqing Education Commission of China (KJQN201801325). Last, but not least, the authors wish to express their gratitude to the anonymous referee whose supporting recommendations provided much help in improving the original manuscript.

\bibliographystyle{mnras}