\section{Introduction}

Traditionally, the determination of the inclination angle, $i$, of a star requires knowledge of the rotation period of the star, \mbox{$P_{\mathrm{rot}}$}\xspace, along with the stellar radius, \mbox{$R_{\ast}$}\xspace, and projected equatorial rotational velocity, \mbox{$\varv\sin i$}\xspace (e.g. \citealt{1994AJ....107..306H, 2007AJ....133.1828W,2010ApJ...719..602S, 2012ApJ...756...66H})
\begin{equation}
    i = \arcsin\left[\frac{\varv \sin i}{\left(2\pi R_{\ast}/P_{\mathrm{rot}}\right)}\right].
\label{eqn:trad_i}
\end{equation}
This technique has been applied to both main-sequence and red-giant stars where such measurements have been available and can offer constraints on the inclination angle \citep{2010ApJ...719..602S, 2012ApJ...756...66H,2015ApJ...807...82T,2017A&A...605A.111C}. Whilst this is not generally a problem for fast rotators for which \mbox{$P_{\mathrm{rot}}$}\xspace or \mbox{$\varv\sin i$}\xspace measurements can be attained, obtaining a measurement of the rotation period or \mbox{$\varv\sin i$}\xspace for slower rotators can prove to be complicated. A difficulty with this method for calculating the inclination angle is that it can lead to unphysical solutions where $\sin i > 1$, which implies that one (or possibly more) of the parameters have been incorrectly inferred, in the context of Eqn~\ref{eqn:trad_i}. The probabilistic method introduced by \cite{2014ApJ...796...47M} aimed to address this by deriving the $\cos i$ distribution given the \mbox{$\varv\sin i$}\xspace, \mbox{$P_{\mathrm{rot}}$}\xspace and radius \mbox{$R_{\ast}$}\xspace distributions. This subsequent technique has been successfully applied to a number of systems and used in larger ensemble analyses \citep{2014ApJ...796...47M, Tiago}.

Asteroseismology offers a means to measure the inclination angle, through the non-radial oscillation modes (e.g. \citealt{1985ApJ...292..238P, 2003ApJ...589.1009G, 2013ApJ...766..101C, 2018arXiv180507044K}). The advantage of asteroseismology is that it can be applied to stars that are relatively slow rotators for which \mbox{$\varv\sin i$}\xspace or \mbox{$P_{\mathrm{rot}}$}\xspace measurements are difficult to obtain. This provides a great opportunity to measure the inclination angle that does not explicitly require the acquisition of \mbox{$\varv\sin i$}\xspace, \mbox{$P_{\mathrm{rot}}$}\xspace or \mbox{$R_{\ast}$}\xspace. Asteroseismology has been used to extract the inclination angle in a number of different situations, such as transiting exoplanetary systems \citep{2013ApJ...766..101C, 2013Sci...342..331H,2014A&A...570A..54L} and in clusters (e.g. \citealt{2017NatAs...1E..64C}). This is the approach we adopt in this work.

There have been many studies where the inclination angle has been an important component, such as in the characterisation of exoplanet and eclipsing binary systems, as it can help reveal the underlying geometry of the system through constraints on the obliquity, the angle between the orbital axis and the stellar rotation axis (e.g. \citealt{1995AcA....45..725G, 2010ApJ...719..602S, 2012ApJ...756...66H,2014ApJ...783....9H,2014ApJ...796...47M,2015MNRAS.446.2959D,Tiago}). In the case of Kepler-56 \citep{2013Sci...342..331H}, the asteroseismically determined inclination angle showed that the system is misaligned (the plane of the planetary orbits are not perpendicular to the rotation axis of the star). As of now, this is the only known misaligned multi-planet system with an evolved host.

In addition to looking at the inclination angle of individual systems, the analysis of an ensemble of inclination angle (or $\sin i$, where $i$ is the inclination angle) measurements is also extremely valuable. By looking at the distribution at a population level, i.e. the distribution of a large sample of stars \citep{2001AJ....122.2008A}, it is possible to test fundamental assumptions used in astrophysical analyses. One such assumption is that the distribution of stellar inclination angles is random (isotropic) with respect to our line of sight \citep{1950ApJ...111..142C}. This is subsequently used in analyses when simulated data are required. For example, the calculation of simulated \mbox{$\varv\sin i$}\xspace values in \cite{2010ApJ...719..602S} involved this assumption to construct the null hypothesis used to test spin-orbit alignment in exoplanet systems. It is also possible to shed light on cluster formation processes through the ensemble analysis of inclination angles. \cite{2010MNRAS.402.1380J} derived the distribution of $\sin i$ values for stars in the open clusters Pleiades and Alpha Per from \mbox{$P_{\mathrm{rot}}$}\xspace, \mbox{$\varv\sin i$}\xspace and radius values derived from the cluster distances, and found no evidence of spin-orbit alignment in either cluster using a Monte-Carlo modelling technique.

The difficulty of extracting the inclination angle has been highlighted by the disagreement in the asteroseismic analyses by \cite{2017NatAs...1E..64C} and \cite{2018arXiv180708301M} when analysing the \emph{Kepler} \citep{2010Sci...327..977B} clusters NGC 6819 and NGC 6791 \citep{2010ApJ...713L.182S, 2011A&A...530A.100H, 2011ApJ...729L..10B, 2011ApJ...739...13S}. \cite{2017NatAs...1E..64C} reported that both clusters showed evidence of strong alignment in their inclination angle distributions, suggesting that the global angular momentum of the initial gas clouds during the cluster formation process were efficiently transferred to the stars leading to a detectable imprint in the inclination angles. However, \cite{2018arXiv180708301M} also analysed these clusters and found that there was no such evidence for strong alignment and in fact that both were consistent with being randomly distributed. This difference in the derived inclination angle arises due to the possibility of misinterpreting whether a star is an extremely slow rotator or possesses a low inclination angle (see section \ref{sec:inc}), which is currently unresolved in these cases.

In this work we follow the method put forward by \cite{2010ApJ...725.2166H} and adopt a Bayesian approach to infer the underlying inclination angle of a set of stars from asteroseismic estimates.

\section[]{Data}

Long-cadence data from the \emph{Kepler} mission are used (a cadence of 29.4 minutes) with all 4 years of observations from quarters Q1-Q17. Our sample consists of 123 stars taken from the 13,000 red giants in \cite{2013ApJ...765L..41S} with a $\nu_{\mathrm{max}}$ (determined by \citealt{2013ApJ...765L..41S}) in the range 231-274 $\mu$Hz. This \mbox{$\nu_{\mathrm{max}}$}\xspace range was chosen because the probability of a star possessing overlapping modes (due to rapid rotation) is very low, therefore making the identification of the mixed modes (those with a mixed p- and g-character, see e.g. \citealt{2014A&A...572L...5M}) easier. All photometric timeseries were produced using the pipeline developed by \cite{2010ApJ...713L.120J}, and power spectra were obtained using the Lomb-Scargle periodogram (\citealt{1976Ap&SS..39..447L,1982ApJ...263..835S}).

Due to the complexity of red giant oscillation spectra, not all of the 123 red giants in our sample could be used in the analysis. A few were observed to show suppressed $\ell=1$ modes \citep{2014A&A...563A..84G, 2015Sci...350..423F, 2016Natur.529..364S, 2017A&A...598A..62M}, which greatly hinder our ability to extract the inclination angle. There were also some cases where the stars had oscillation modes reflected from across the Nyquist frequency \citep{2014MNRAS.445..946C}. This made the mode identification difficult because there was a mixture of real and reflected oscillation modes obscuring the well established frequency patterns. Finally, there were also a few stars for which the modes could not be successfully disentangled with the current method due to the highly complex spectra.

\section[]{Data Analysis}
\label{sec:analysis}

\subsection{Mode Detection}\label{sec:det}

The first step towards extracting the inclination angle is to detect the $\ell=1$ mixed modes to which we want to fit our model. For the detection of the modes we follow the method given in \cite{2004A&A...428.1039A} applied over each radial order. Rather than applying the detection test to one realisation of the data, we opt to apply the test to a number of rebinned realisations (30 in total, varying from 2 to 60 bins). This amounts to a resolution of $\sim0.03\Delta\nu$ at the high $\nu_{\mathrm{max}}$ values of the stars in our sample. The purpose of including the rebinned realisations is to try and accentuate the features of narrow modes as opposed to noise spikes. A false alarm probability of 10\% is used over a window in frequency of $\Delta\nu/2$, under the assumption of having $\Delta\nu/\left(2N\nu_{\mathrm{bw}}\right)$ independent bins (where $N$ is the number of bins binned over and $\nu_{\mathrm{bw}}$ is the frequency resolution of the rebinned realisation). The detection test is performed for each rebinned realisation and the frequencies of each realisation that exceed the detection threshold are kept. These frequencies are then clustered using the mean shift clustering algorithm \citep{Comaniciu:2002:MSR:513073.513076}\footnote{As implemented in the \texttt{python} package \texttt{scikit-learn} \citep{scikit-learn}.}. We set the minimum number of points needed to constitute a cluster to 5, ensuring that single noise spikes in a few realisations are not kept. The value of 5 was chosen because we empirical found that it was a good trade-off between detecting narrow modes with lower signal-to-noise ratios and detecting spikes due to noise. The identified cluster centres are taken as initial guesses for the mode frequencies in the fitting process (see section~\ref{sec:inc}). In addition, initial guesses for the rotational splitting and inclination angle can also be approximated from the cluster centres. It is important to note that we only search for modes between the $\ell=0$ and $\ell=2$ mode of the same order to avoid ambiguity and issues in the fitting procedure. This is why, in Fig.~\ref{fig:detection}, the mixed mode at just above 245 $\mu$Hz is not used. An example of the mixed modes detected by this algorithm is shown in Fig.~\ref{fig:detection}.

\subsection{Mode Identification}\label{sec:ID}

The detected peaks are of no use by themselves and in order to fit our model to the data we need to know which peaks belong to which mixed mode and whether they are rotationally split. The initial identification of the radial modes is performed using the universal pattern described by \cite{2011A&A...525L...9M}, which is dependent only upon the large frequency separation, $\Delta\nu$.

We follow the formulation for $\ell=1$ mixed modes frequencies given in \cite{2012A&A...540A.143M}

\begin{figure*}
\centering
\includegraphics[width=\textwidth]{Figure1.pdf}
\caption{The top panel shows a log-log plot of the power spectrum for KIC 5553307 in black and a smoothed spectrum ($1\;\mu$Hz boxcar) in red. The centre and bottom panels show an example radial order taken from the same star. The centre panel shows the modes that have been detected following the method given in section~\ref{sec:det}. The cyan and magenta triangles denote the possible positions of the detected rotationally split components ($m=-1$ and $m=+1$ respectively). Only the $m=\pm1$ components are seen due to the inclination angle of the star being close to 90\mbox{$^{\circ}$}\xspace. The red dashed line in the centre panel shows the detection threshold assuming a false-alarm probability of 10\%. The bottom panel shows the modes chosen once the additional constraint of the mode identification (see Section~\ref{sec:ID}) has been applied. For example, the peaks at $\sim$230$\mu$Hz, whilst detected, were not used in the subsequent analysis as they do not agree with the assumed mixed mode pattern.}
\label{fig:detection}
\end{figure*}

\begin{equation}
\nu  = \nu_{n_{p},\ell=1} +  \frac{\Delta\nu}{\pi}\arctan\left[q\tan\pi\left(\frac{1}{\Delta\Pi_{1}\nu} - \varepsilon_{g}\right)\right],
\label{eqn:asymp}
\end{equation}
where $\nu_{n_{p},\ell=1}$ is the nominal p-mode frequency, $\Delta\nu$ is the large frequency separation, $q$ is the coupling between the p- and g-modes, $\Delta\Pi_{1}$ is the $\ell=1$ period spacing and $\varepsilon_{g}$ is a phase term. It is assumed that $\varepsilon_{g}=0$, in accordance with \cite{2011Natur.471..608B} and \cite{2012A&A...540A.143M}. This assumption has since been shown to be questionable by \cite{2018A&A...610A..80H}. Given that the parameters were used only to guide the eye this is nevertheless justifiable, since the procedure was more than ample to allow robust identification of the modes.

Eq.~\ref{eqn:asymp} is combined with the formulation for the expected rotational splittings of $\ell=1$ mixed modes as described in \cite{2013A&A...549A..75G}

\begin{equation}
\frac{\nu_{\mathrm{s}}}{\nu_{\mathrm{s}, \mathrm{max}}} = \zeta\left(1-2\frac{\left\langle\Omega\right\rangle_{\mathrm{env}}}{\left\langle\Omega\right\rangle_{\mathrm{core}}}\right) + 2 \frac{\left\langle\Omega\right\rangle_{\mathrm{env}}}{\left\langle\Omega\right\rangle_{\mathrm{core}}},
\label{eqn:splittings}
\end{equation}
where $\nu_{\mathrm{s}}$ is the rotational splitting, $\nu_{\mathrm{s}, \mathrm{max}}$ is the maximum splitting, $\zeta$ is the ratio of the mode inertia in the g-mode cavity to that of the entire cavity \citep{2015A&A...580A..96D}, $\left\langle\Omega\right\rangle_{\mathrm{core}}$ is the angular rotational velocity averaged over the core regions and $\left\langle\Omega\right\rangle_{\mathrm{env}}$  is the angular rotational velocity averaged over the envelope. We make the assumption that the contribution from the envelope is very small, i.e. that the envelope is rotating very slowly. In the limit $\left\langle\Omega\right\rangle_{\mathrm{env}} \ll \left\langle\Omega\right\rangle_{\mathrm{core}}$ the ratio $\left\langle\Omega\right\rangle_{\mathrm{env}}/\left\langle\Omega\right\rangle_{\mathrm{core}}$ can be neglected, reducing Eq.~\ref{eqn:splittings} to

\begin{equation}
\nu_{\mathrm{s}} = \zeta\nu_{\mathrm{s}, \mathrm{max}}.
\end{equation}

The individual parameters contained in equations \ref{eqn:asymp} and \ref{eqn:splittings} were manually adjusted to produce a good fit by eye. The values derived are approximate and are only used in the generation of the artificial data; otherwise they are not required in the rest of the inference. In other words, we construct a pattern of the approximate $\ell=1$ mixed mode central frequencies and rotational splittings to identify the mixed modes. The mixed mode is selected if we observe the prerequisite number of peaks: one for angles close to zero, two for those close to 90\mbox{$^{\circ}$}\xspace and 3 for intermediate angles. These provide initial guesses for the central frequency and the rotational splitting of the mixed modes, which are then subsequently fed into the fitting procedure. An example of the final modes selected is given in Fig.~\ref{fig:detection}. It might be expected that by only fitting a subset of our data (due to not detected all components of a mode) that we could introduce a selection effect. However, an underlying assumption that is made during the subsequent analysis is that each mode possesses the same underlying  inclination angle. As a result, this selection effect will only impact the uncertainty of the determination of the inclination angle rather than the underlying value itself.

\subsection{Peakbagging Model}\label{sec:inc}

The determination of the inclination angle uses the formulation described in \cite{2003ApJ...589.1009G}. This relies on the assumption that there is equipartition of energy between mode components of differing azimuthal order in any given multiplet. For $\ell=1$ modes the inclination angle can be derived from the amplitudes of the $|m|=\ell$ and $m=0$ components.

The model that is fitted to each mixed mode is given by (e.g. \citealt{refId0}):
\begin{equation}
\mathcal{M}(\nu; \mathbf{\theta}) = \sum^{\ell}_{m=-\ell} \frac{\mathcal{E}_{1m}(i)H}{1 + \frac{4}{\Gamma^{2}}(\nu-\nu_{0}-m\nu_{s})^{2}}+ B.
\label{model}
\end{equation}
where $B$ is the background (the very narrow region in frequency in which the mixed mode is being fitted results in the use of a flat background being valid), $\nu_{0}$ is the central frequency of the mixed mode (i.e. the frequency of the $m=0$ component), $\nu_{s}$ is the rotational splitting and the summation is over each $m$ component which runs from $-\ell$ to $+\ell$. The height $H$ can be parametrised in terms of the amplitude of the mode, $A$, and the linewidth, $\Gamma$, to alleviate the unwanted impact of the strong anti-correlation between $H$ and $\Gamma$ \citep{2003A&A...398..305C, 2006MNRAS.371..935F}
\begin{equation}
    H = \frac{2A^{2}T}{\pi\Gamma T + 2},
\end{equation}
where $T$ is the length of the observations. This expression is used in order to account for the change in the mode profile when the linewidth tends towards the unresolved regime, i.e. $\Gamma \lesssim 2\nu_{\mathrm{bw}}/\pi$.

The final parameter to introduce is the visibility factor $\mathcal{E}$ \citep{2003ApJ...589.1009G} which is given by
\begin{equation}
\mathcal{E}_{\ell m}(i) = \frac{(\ell - |m|)!}{(\ell+|m|)!}P^{|m|}_{\ell}(\cos i)^{2},
\label{vis}
\end{equation}
subject to the constraint
\begin{equation}
\sum_{m}\mathcal{E}_{\ell m}(i) = 1,
\end{equation}
where $P^{|m|}_{\ell}(\cos i)^{2}$ are associated Legendre polynomials. For the case of $\ell=1$ modes, Eq~\ref{vis} gives
\begin{align}
\mathcal{E}_{1,0}(i) &= \cos^{2}(i),\\
\mathcal{E}_{1,\pm1}(i) &= \frac{1}{2}\sin^{2}(i).
\label{eqn: inc_leg}
\end{align}

The mode visibility therefore modulates the amplitudes of the modes in a way that is dependent upon the inclination angle of the star. If we were to observe a star at 90\mbox{$^{\circ}$}\xspace (i.e. equator-on) with respect to our line of sight then in the power spectrum we would only see the outer components of the $\ell=1$ mixed modes (i.e. $m=\pm1$). Whereas if we were to observe a star pole-on, at close to 0\mbox{$^{\circ}$}\xspace with respect to our line of sight, then in the power spectrum we would only observe the central ($m=0$) component. Intermediate angles would fall in between the two extremes in terms of relative amplitudes. An example of how the inclination angle modulates the component amplitudes of an $\ell=1$ mode is shown in Fig.~\ref{fig: splittings_nice}.

\begin{figure*}
\centering
\subfloat{\includegraphics[width=0.8\textwidth]{Figure2.pdf}}
\caption{Example $\ell=1$ modes observed at different inclination angles. The $m=\pm1$ components are shown in orange and the central $m=0$ component is shown in blue. The angle of inclination used in simulating the mode is labelled in each panel. Furthermore, a mode linewidth of $0.2\mu$Hz, a rotational splitting of $0.4\mu$Hz and an observing length of 4 years were used in all cases.}
\label{fig: splittings_nice}
\end{figure*}

\subsection{The Fitting Process}\label{sec:fit}

Our interest lies in fitting individual rotationally split $\ell=1$ mixed modes rather than performing a global fit. There are many mixed modes per order and the high quality \emph{Kepler} data enable us to extract many individual measures of the inclination angle.

In this work we adopt a Bayesian framework for fitting the mixed modes (see e.g \citealt{refId0,2014A&A...570A..41K,2016MNRAS.456.2183D}) since we want to obtain posterior probability distributions for the inclination angle to be used in the subsequent analysis. We use Markov Chain Monte Carlo (MCMC) to sample the parameter space, making use of \texttt{ptemcee}\footnote{\texttt{ptemcee} is a parallel tempering extension to \texttt{emcee} that uses dynamic temperature selection.} \citep{2016MNRAS.455.1919V, 2013PASP..125..306F}. In subsequent sections we will adopt this procedure as well.

The log-likelihood function used in the model fitting assumes the noise properties of the power spectrum follows a $\chi^{2}_{2}$ distribution and is given by \citep{1986ASIC..169..105D, 1990ApJ...364..699A}
\begin{equation}
\ln L = -\sum_{\nu}\left[\ln M_{\nu}(\mathbf{\theta}) + \frac{P_{\nu}}{M_{\nu}(\mathbf{\theta})}\right],
\end{equation}
where the summation is made over each frequency bin, $M_{\nu}(\mathbf{\theta})$ is the model evaluated at a given frequency for a set of parameters $\mathbf{\theta}$, and $P_{\nu}$ is the power at a given frequency.

The prior distributions were taken to be uniform for all parameters apart from the central frequency of the mode and the inclination angle. A Gaussian prior was placed on the central frequency with a mean according to the approximate mode frequency extracted from the clustering (see section~\ref{sec:det}) and a standard deviation of 0.2 $\mu$Hz. This is a weakly informative prior, the standard deviation of which is taken to be a value that is narrow enough that the central frequency is fitted to the central component of the mixed mode and wide enough to allow for the fact that the initial guess may not exactly coincide with the underlying value. Finally an isotropic prior, $p(i) \propto \sin i$, was placed on the inclination angle which comes from the assumption that stars are oriented randomly with respect to the observer. The isotropic prior is also uninformative and preferable to a uniform prior, as shown in more detail in Appendix~\ref{sec: why_uninform}. The prior distributions used are given in Table~\ref{tab:priors}.

\begin{table}
	\centering
	\caption{Model prior distributions for the peak-bagging analysis. $\mathcal{N}$(mean,standard deviation) indicates Gaussian priors, and $\mathcal{U}$(lower bound, upper bound) indicates uniform priors.}
	\label{tab:priors}
	\begin{tabular}{ll}
	\hline
	Parameter & Prior\\
	\hline
	$A$ & $\mathcal{U}(0,50)$ (ppm)\\
	$\nu_{0}$ &$\mathcal{N}(\nu_{0,\mathrm{initial}},0.2)$ ($\mu$Hz)\\
	$\Gamma$ & $\mathcal{U}(0,5)$ ($\mu$Hz)\\
	$\nu_{\mathrm{s}}$ & $\mathcal{U}(0, 0.7)$ ($\mu$Hz)\\
	$i$ & $\sin i$ (degrees)\\
	$B$ & $\mathcal{U}(0, 200)$ (ppm$^{2}\mu$Hz$^{-1}$)\\
	\hline
	\end{tabular}
\end{table}

\begin{figure}
\centering
\subfloat{\includegraphics[width=0.45\textwidth]{Figure3.pdf}}
\caption{An example of two Lorentzian fits made to modes in KIC 5553307 are shown in the left-hand side and their respective inclination angle posterior probability distributions are shown on the right-hand side. The top panel shows a more g-dominated $\ell=1$ mixed mode at $235.5\;\mu$Hz, and the resolved p-dominated $\ell=1$ mixed mode at $238.3\;\mu$Hz is shown in the bottom panel. On the left-hand side, the power spectrum is shown in black with the best fit model overplotted in red. The respective $m$ components are also labelled as well as a demonstration of the rotational splitting $\nu_{s}$. The black dotted line denotes the estimated position of the $m=0$ component which is not visible due to the close to 90\mbox{$^{\circ}$}\xspace angle of inclination of the star.}
\label{fig: splittings}
\end{figure}

This approach yields posterior distributions of the inclination angle for individual modes rather than point estimates, which will be useful in the following sections. An example of fits made to both p- and g-dominated modes is shown in Fig.~\ref{fig: splittings}.

\section[]{Extracting the inclination angle}\label{sec: method1}

\subsection{Hierarchical Method}

We adopt an approach in this work that explicitly models the underlying distribution whilst accounting for the initial prior distribution. This is performed following the hierarchical Bayesian method described in \cite{2010ApJ...725.2166H} and again using MCMC (see section~\ref{sec:fit}) to sample the parameter space. The explicit modelling of the distribution should enable us to better approximate its underlying shape, especially in the case of truncation at low and high inclination angles. The method is hierarchical because we are allowing the prior distribution to be modelled. We are essentially trying to find a model that better describes the data than the original uninformative prior. This enables us to infer the underlying inclination angle from the model given the posterior samples that we gathered in the previous section.

In order to infer the underlying angle, we turn to the marginalised likelihood given in \cite{2010ApJ...725.2166H} (see Appendix~\ref{sec:like_deriv} for a full derivation)

\begin{equation}
\mathcal{L_{\alpha}} \approx \prod^{N}_{n=1}\frac{1}{K}\sum^{K}_{k=1}\frac{f_{\alpha}(i_{nk})}{p_{0}(i_{nk})}.
\label{eqn: new_like}
\end{equation}
where $f_{\alpha}(i_{nk})$ is the model we want to fit with parameters $\alpha$, evaluated at the $k$th sample of the posterior of the $n$th mode. The parameter $p_{0}(i_{nk})$ is the uninformative prior (used in the original fitting process, as denoted by the subscript 0) evaluated for the same sample.

In addition, Eq.~\ref{eqn: new_like} can be interpreted as follows: inside the sum is the ratio of the new prior to the initial (uninformative prior), taken over each mode. The prior used in the initial fitting process should be uninformative and be valid over the entire range occupied by the new prior, otherwise the inferred parameters can be very uncertain. Eq.~\ref{eqn: new_like} is a marginalised likelihood because we have integrated out all of the parameters of the original fit. We have therefore assumed that the distribution over the parameters from the original fit is separable, which has enabled us to formulate the likelihood function in terms of just the inclination angle. By assuming separability, we have also assumed that the only distribution that needs to be adjusted from the assumed prior to match the population distribution is the inclination angle.

We adopt the hierarchical method because it allows us to combine the multiple observations of inclination angle (one per mixed mode) in a principled way, accounting for the uncertainty in each observation and also possible systematics in the amplitudes and rotational splitting of each mode. In order to model the systematics, we do not assume that the inclination parameters measured by the fit to each mode are identical, we rather assume that they come from a distribution with a central peak and scatter; see Section~\ref{sec:model}. By fitting for both the central value and the scatter (i.e. the location and concentration parameters defined in Section~\ref{sec:model}), we can obtain an estimate of the inclination angle and the degree of systematic scatter in the mode-to-mode values simultaneously, all while accounting for the uncertainties in each mode's estimate of the inclination angle.

\subsection[]{Inclination angle model of the hierarchical method}\label{sec:model}

The model we adopt for $f_{\alpha}(i)$ is a slightly modified version of the Fisher distribution (as chosen in \citealt{2009ApJ...696.1230F}) that is equivalent to a zero-mean Normal distribution on the sky. Whilst the original formulation of the Fisher distribution is a valid model if we assume our angle of inclination distribution is (or is very close to) isotropic, it is unsuitable for the highly localised inclination angle distributions of individual stars.

To add flexibility to the model, a location parameter is added that enables the peak of the distribution to be shifted in the region $\mu\in[0, \pi/2]$. This results in the model being able to adapt to both a sharply peaked distribution (in the case of the individual stars) and a much wider population distribution. It offers the ability to model an isotropic distribution and one localised towards any angle between 0 and 90 degrees (i.e. towards anisotropy). The location parameter will provide information on the underlying inclination angle of the star and so this is our parameter of interest. The updated form is given by
\begin{equation}
f(i|\mu,\kappa) = \exp\left[\kappa\cos\left(i - \mu\right) \right]\sin i,
\label{eqn: model}
\end{equation}
where $\kappa$ is the concentration parameter and the additional parameter $\mu$ is the aforementioned location parameter. We have used $f$ to explicitly state that the function is not a probability distribution, and will proceed to use $p$ otherwise. In its current state the distribution in Eq.~\ref{eqn: model} is unnormalised.  The normalisation constant can be derived analytically (as shown in Appendix~\ref{sec: norm_deriv}) to give the probability distribution
\begin{equation}
p(i|\mu,\kappa) = \left\{I_{0}(\kappa) + 2\varphi(\kappa,\mu)\right\}^{-1}
\exp\left[\kappa(i-\mu)\right]\sin i,
\end{equation}
where $I_{0}(\kappa)$ is a modified Bessel function of the first kind of order zero, and $\varphi(\kappa,\mu)$ is a function of the Fisher distribution parameters used in the normalisation of the distribution (for more information see Appendix~\ref{sec: norm_deriv}).

Examples of this modified Fisher distribution, which will become our model in the hierarchical method, are shown in the top panel of Fig.~\ref{fig: new_fisher}. The distribution can represent isotropy when $\kappa=0$, and the addition of the location parameter also allows for deviations from isotropy to be accounted for. The average measured inclination angle under the assumption of isotropy is $\sim57.2$ degrees (1 radian) and so the expected value of $\mu$ should also take the same value under isotropy\footnote{The location parameter $\mu$ will be fitted in radians and so will be referred to and displayed as such, otherwise we will refer to the angles in degrees.}. Therefore, subtle deviations or possible biases could be inferred from this parameter and not just from $\kappa$. In addition, the standard deviation of the distribution as a function of $\kappa$ is given in the last panel of Fig.~\ref{fig: new_fisher}. This helps show how the variance of the distribution decreases with increasing $\kappa$ leading to a much more localised and sharply peaked distribution.

An important consideration in the hierarchical analysis is the support (i.e. region of validity in parameter space) of both the initial prior and the new prior model in Equation~\ref{eqn: new_like}. It is advisable to avoid priors that tend to zero at the edges of the parameter space since this can cause problems in the likelihood function (Eq.~\ref{eqn: new_like}), due to the division. This occurs if the model is fitted in inclination angle space due to the isotropic prior tending to zero with the inclination angle. However, this effect can be mitigated by instead performing the inference in $\cos i$ whereby the isotropic prior becomes uniform (see Appendix~\ref{sec: why_uninform}). Therefore the modified Fisher distribution as a function of $i$ (as given above) is not adequate in the modelling and so must be transformed such that we obtain the distribution in $\cos i$, given that $i$ is Fisher distributed with some concentration parameter $\kappa$ and location parameter $\mu$.

Following the method given in \cite{2014ApJ...796...47M} for transforming the original Fisher Distribution, we can transform this updated model (as shown in Appendix~\ref{sec: model}) to derive the following probability distribution

\begin{equation}
p(y | \mu,\kappa) = \left\{I_{0}(\kappa) + 2\varphi(\kappa,\mu)\right\}^{-1}
\exp\left(\kappa y\cos\mu + \kappa\sqrt{1-y^{2}}\sin\mu\right),
\label{new_fisher}
\end{equation}
where $y = \cos i$.

\begin{figure}
\centering
\subfloat[]{\includegraphics[width=0.45\textwidth]{Figure4a.pdf}}
\vfill
\subfloat[]{\includegraphics[width=0.45\textwidth]{Figure4b.pdf}}
\vfill
\subfloat[]{\includegraphics[width=0.45\textwidth]{Figure4c.pdf}}

\caption{A demonstration of the model used in the hierarchical modelling as a function of $i$ and $\cos i$. The top panel shows the modified Fisher distribution as a function of $i$ for a variety of different $\kappa$ values. The solid lines denote distributions computed with $\mu=\pi/2$ and the dashed lines show distributions computed with $\mu=0$. The middle panel shows the distributions in the top panel as a function of $\cos i$ using equation \ref{new_fisher}. The bottom panel shows how the parameter $\kappa$ affects the standard deviation of the distribution, with larger $\kappa$ leading to a narrower distribution.}
\label{fig: new_fisher}
\end{figure}

An isotropic prior was imposed on the location parameter $\mu$ in order to reflect the isotropic prior distribution on the angle during the fitting process. The application of an uninformative prior on $\mu$ is important, as we shall see in the next section. When inferring the inclination angle distribution for individual stars, the distributions are expected to be much narrower than the population distribution (since we are searching for localised distributions). As a result, large values of $\kappa$ are to be expected and the prior adopted by \cite{2014ApJ...796...47M},  $p(\kappa) = (1 + \kappa^{2})^{-3/4}$, is too aggressive and prevents large $\kappa$ values from being explored. Therefore, we instead adopt a Half-Cauchy prior with a width, $\gamma=50$, on $\kappa$ (e.g. \citealt{2009arXiv0901.4011G}):

\begin{equation}
p(\kappa | \gamma) = \frac{1}{\pi\gamma}\left(\frac{\gamma^{2}}{\kappa^{2} + \gamma^{2}}\right)
\;\text{ for}\;\kappa>0.
\end{equation}
The extended tail of the Half-Cauchy distribution enables much larger $\kappa$ values to be explored if necessary. The width of the Half-Cauchy prior was tested using values from 10 to 100 and in all cases consistent values were returned indicating that the choice of $\gamma$ is not adversely affecting the inferred parameters.

\section{Application to data}

\subsection{Artificial data}

The method was first tested on artificial red giant power spectra with properties drawn from the same parameter space as the real data and with inclination angles drawn from a known distribution (i.e. isotropic). This testing shows the expected efficiency when applied to real data, such that the accuracy and precision of the method can be assessed. Where available, the same parameters are used in the construction of these power spectra (e.g. estimated period spacings, coupling factors etc.). The construction of the artificial datasets is discussed in Appendix~\ref{sec: artif_spec}\footnote{The code used to generate these power spectra can be found at \url{https://github.com/jsk389/Artificial-red-giant-power-spectra}.}.

The mixed modes in the power spectrum are identified and fitted in the manner described in Section~\ref{sec:analysis}, resulting in inclination angle posterior distributions of every fitted mode for each star in the sample. A total of 77 of the 90 stars in the simulated sample rotationally split modes were fitted properly. This is because the rotational splitting values are drawn from the estimated rotational splittings extracted whilst performing the mode identification on the real stars. Since there are stars for which the rotational splitting could not be observed (either because of an inclination angle of zero or being a very slow rotator) their rotational splitting estimates were set to zero. As a result there was a chance of drawing a non-zero inclination angle from the input distribution and drawing the rotational splitting parameter from a star with no observed rotational splitting. This would result in an inferred inclination angle inconsistent with the input angle.

Examples of the application of the hierarchical analysis to a few stars in the artificial sample are shown in Fig.~\ref{fig: artif_example_1}. A Monte-Carlo analysis, performed by computing the weighted mean, is also shown for comparison. The weighted mean is computed by drawing one value from the inclination angle distribution of each mode and calculating the mean of the drawn values weighted by the inverse of the variance of the distribution it was drawn from. This is then repeated 1000 times to build up a distribution of the weighted mean. In the top panel, the extended tails towards 90\mbox{$^{\circ}$}\xspace bias the Monte-Carlo estimate away from the input angle, although this is only slight for an input angle close to 50\mbox{$^{\circ}$}\xspace (where such tails should not occur often). This highlights the benefit of using the hierarchical inference, since the effect of the prior, which, as discussed earlier, is a major contributor to these tails, is greatly reduced and the posterior PDF of $\mu$ is nicely centred about the input angle. A more extreme case is shown in the bottom panel of Fig.~\ref{fig: artif_example_1} where the input angle is close to 90\mbox{$^{\circ}$}\xspace and the posterior PDFs of the inclination angle show more structure that is indicative of a lower signal-to-noise fit. Even in such a case, it is reasonably clear as to approximately where the input angle should lie. The Monte-Carlo estimate has been biased by the truncation of the data at 90\mbox{$^{\circ}$}\xspace.

The results of the hierarchical method applied to the full artificial dataset is shown in Fig.~\ref{fig: artif_comparison}. In the vast majority of cases, the inferred inclination angle (i.e. the posterior distribution of the location parameter) is consistent with the input angle from the posterior estimates. This is promising and shows that we can recover the input angle and in fact all but one of the stars in the sample lies within the 95\% highest posterior density (HPD) credible interval of the underlying value which shows excellent agreement. The results from the hierarchical analysis verify our interpretation of the updated model that the location parameter reflects the underlying angle we wish to infer.

The hierarchical model is not perfect however, and it is important to explain the apparent structure present for inclination angles above 80\mbox{$^{\circ}$}\xspace and at low angles. Both can be explained by the fact that the parameter $\kappa$ becomes much more important close to the edge of the region of support, i.e. close to 0\mbox{$^{\circ}$}\xspace or 90\mbox{$^{\circ}$}\xspace. For large $\kappa$ values ($\kappa > 50$, i.e. a narrow distribution), there is degeneracy between $\kappa$ and $\mu$ when $\mu\rightarrow\pi/2$ and so even a large change in $\kappa$ will not result in a significant change in the shape of the distribution. This uncertainty will propagate through into the posterior PDF of $\mu$. In the 90\mbox{$^{\circ}$}\xspace case this is evident, whereby the posterior peaks very close to 90\mbox{$^{\circ}$}\xspace and we cannot say any more precisely as to where the mean value should lie. This idea will be expanded upon in section~\ref{sec:discussion}.

For comparison we also include Monte-Carlo estimates of the weighted mean of the inclination angle posterior distributions, which are shown in Fig.~\ref{fig: artif_comparison}. It is clear that the Monte-Carlo estimate performs poorly close to 0\mbox{$^{\circ}$}\xspace and 90\mbox{$^{\circ}$}\xspace as expected due to the effects of the prior and truncation respectively. In the approximate region of 30\mbox{$^{\circ}$}\xspace to 70\mbox{$^{\circ}$}\xspace, both the Monte-Carlo estimate and the hierarchical analysis give very similar results. The poorer performance of the Monte-Carlo method is evident by the fact that 22 stars are not within the 95\% HPD credible interval, which occurs at very low ($i<20^{\circ}$) and high ($i > 80^{\circ}$) angles.

The reasons for the Monte-Carlo estimate failing to accurately calculate the underlying inclination angle across all of parameter space can be summarised as follows. Firstly, when the inclination angle is difficult to constrain, the prior can have a non-negligible contribution to the posterior distribution. This is especially apparent because our prior on the inclination angle is not uniform and a large amount of probability mass is contained close to 90 degrees. As a result the posterior distributions in angle of lower signal to noise modes can have an extended tail towards higher angles. The presence of such a feature would cause a bias in both the mean and variance of the Monte-Carlo estimate. Correcting the effect of the prior \emph{a posteriori} should be avoided and it would be much better to account for this at the same time as approximating the underlying angle. Secondly, the Monte-Carlo estimate fails as a result of the way we measure the inclination angle of the star. Due to the nature of the inclination angle and its symmetries (e.g. \citealt{2003ApJ...589.1009G, 2017NatAs...1E..64C}), we cannot distinguish between an angle $i$ or $(180 - i)^{\circ}$ and so our inclination angles are restricted to lie with the range $i\in [0, 90]^{\circ}$. This introduces the problem of truncation near $i=90^{\circ}$ that will cause the Monte-Carlo estimate to bias low away from $i=90^{\circ}$ due to where the mean lies for the truncated distribution. Additionally this method will bias high when the angle of inclination is close to $i=0^{\circ}$, due to the extended tail seen in the posterior probability density functions (PDFs) of low inclination angle stars. This highlights a need to account for the presence of the prior, which in turn is, to a large extent, responsible for those extended tails. Examples of this can be seen in Fig.~\ref{fig: artif_example_1}.

\begin{figure}
\centering
\subfloat[]{\includegraphics[width=0.45\textwidth]{Figure5a.pdf}}
\vfill
\subfloat[]{\includegraphics[width=0.45\textwidth]{Figure5b.pdf}}
\caption{Examples of the fitting procedure applied to two stars from the artificial sample, the top panel shows a star with an input angle (denoted by the red dashed line) of $\sim50^{\circ}$ and the bottom panel shows a star with an input angle close to 90\mbox{$^{\circ}$}\xspace. The individual overlapping posteriors for each mode are plotted in blue. The mean inclination angle inferred from the Monte-Carlo estimate is shown in orange and the posterior PDF of the location parameter $\mu$ is shown in black.}
\label{fig: artif_example_1}
\end{figure}

\begin{figure}
\centering
\includegraphics[width=0.45\textwidth]{Figure6a.pdf}
\caption{The inferred inclination angle shown against the input inclination angle for the artificial data generated with an underlying isotropic distribution. The inferred angles from the hierarchical analyses are shown in black, and those from the Monte-Carlo estimate are shown in orange. In both cases the error bars are taken to be the 68.3\% highest posterior density (HPD) interval. The residuals of the data about the 1:1 line (the red dotted line) are shown at the bottom.}
\label{fig: artif_comparison}
\end{figure}

\subsection{Application to real data}

We applied the same method described above for mode identification, detection and fitting to the \emph{Kepler} red giants in our sample, resulting in posterior PDFs for the inclination angle of each star. The sample dropped from 90 to 89 stars as for one star only a single mode could be fitted and there is no means to check the consistency of the fit. A table of the inferred inclination angles are given in Table~\ref{tab:real_incs}.

\section[]{Population Distribution}

Now that we have obtained inclination angles for individual stars, we can also go one step further and assess whether or not this population of stars is isotropically (randomly) distributed in angle. To do so we opt for a straightforward approach based around the Kolmogorov-Smirnov (K-S) test (see e.g. \citealt{doi:10.1080/01621459.1951.10500769}).

First of all, we calculate the empirical cumulative distribution functions (ECDFs) of the artificial and real data by drawing an inclination angle from each of the inclination angle distribution for the stars we obtained a measurement for, constructing the ECDFs and then repeating this 1000 times, as shown in Fig~\ref{fig:cdfs}. Inspection of the ECDFs reveals the angles at which excesses lie compared to an underlying isotropic distribution. The gradient of the ECDF provides information about these excesses and where they occur, for example, at angles close to 0\mbox{$^{\circ}$}\xspace the ECDF of both datasets rises faster than the isotropic distribution reference indicating an excess and the subsequent plateau shows a lack of stars with respect to the isotropic distribution. For the intermediate angles there is good consistency between the distributions, indicating that the isotropic distribution represents the data well in this region. For high angles the ECDF of the datasets rise much faster than the isotropic distribution, again indicating an excess at those high angles. Note that this point is approximately where the largest deviation is found between the ECDFs and the isotropic distribution.

For the K-S test, the isotropic distribution is taken as our reference distribution and what we wish to understand is what the expected distribution of stars would look like given the number of stars in our sample. To build our reference, we randomly draw an inclination angle (of the 89 stars in the case of the real data and 77 stars in the case of the artificial data) from an isotropic distribution and then calculate the empirical cumulative distribution function (ECDF) of those draws. We then compute the maximum deviation $D$ between the ECDF and the CDF of the isotropic distribution $P(i)\propto 1-\cos i$. This is repeated $1\times10^{5}$ times to build up a reference distribution of the maximum deviation expected for an isotropic distribution of our sample size, as shown in Fig.~\ref{fig:d_stat}.

The distribution of maximum deviation gives us the ability to assess where both the real data and artificial data lie in terms of consistency with an underlying isotropic distribution for the given sample size, this is highlighted in Fig~\ref{fig:cdfs} by the red lines. In the case of the artificial data the initial inclination angles were drawn from an isotropic distribution and so, as expected, the value of $D$ is consistent with this, lying at the 79th percentile. For the real data we do not know the underlying distribution, and the distribution of $D$ lies at the 99th percentile.

The discrepancy seen in the maximum deviation is unexpected. Given the way we have selected the stars in our sample there is no anticipated reason for the population distribution not to be consistent with isotropy. Before we can draw conclusions about the underlying population distribution, however, we must also consider that there are observational and population level biases that could well explain this discrepancy, such as the artificial simulations lacking physics that is present in the real data. These biases need only be small, and in individual cases can be less than the uncertainty on the inclination angle, but they are compounded when looking at the population as a whole. This motivates the need for a selection function.

The fundamental quantity required to estimate the effect of the selection function on the population inference is the average detectability of a member of the population $\langle P_{\mathrm{det}}(i)\rangle$ \citep{2019MNRAS.486.1086M}
\begin{equation}
	\langle P_{\mathrm{det}}(i)\rangle = \int P_{\mathrm{det}}(i)p(i|\lambda)\mathrm{d}i,
\end{equation}
where subscript $\mathrm{det}$ denotes stars that have a detected inclination angle and $\lambda$ are the parameters of the population distribution (e.g. a modified Fisher distribution as used in this work, whose parameters can be inferred using an extension of the hierarchical inference to the population of stars). Due to the ``manual'' nature of the pipeline we use to estimate the inclination angle, estimating this quantity is beyond the scope of this paper; if the pipeline could be automated (e.g. \citealt{2018MNRAS.476.1470G}), then \cite{2019RNAAS...3e..66F} provides a simple way to estimate this quantity via a Monte-Carlo method with synthetic populations of stars with varying inclination angle and noise properties.

\begin{figure}
\centering
\subfloat[]{\includegraphics[width=0.45\textwidth]{Figure7a.pdf}}
\vfill
\subfloat[]{\includegraphics[width=0.45\textwidth]{Figure7b.pdf}}
\caption{Empirical cumulative distribution functions (ECDFs) for the artificial data (top-panel) and real data (bottom panel) shown for 1000 draws (as described in text). The percentiles of the distribution are given along with the isotropic distribution shown by the black dashed line. The position of the largest deviation in each case is given by the red line.}
\label{fig:cdfs}
\end{figure}

\begin{figure}
\centering
\centering
\subfloat[]{\includegraphics[width=0.45\textwidth]{Figure8a.pdf}}
\vfill
\subfloat[]{\includegraphics[width=0.45\textwidth]{Figure8b.pdf}}
\caption{The distribution of the D-statistic for our sample size and $1\times10^{5}$ random draws. The top panel shows the distribution for the artificial data and the bottom panel shows the distribution for the real data, each time the D-statistic computed for each dataset is given by the dashed line.}
\label{fig:d_stat}
\end{figure}

\section[]{Discussion}\label{sec:discussion}

In order for the hierarchical inference to work in the way presented here, we have made the assumption that all priors are separable and there are no dependencies of inclination angle on other parameters. This is weakly violated as the inclination angle is inferred through the relative mode heights. Therefore, if we were to make this method truly hierarchical, we would have to resample the amplitudes as well. This is, however, not necessary since both the artificial data and real data occupy the same region of parameter space in terms of signal-to-noise ratio and reduced splitting (the ratio of mode linewidth to the rotational splitting) as shown in Fig.~\ref{fig: red_split_comp}. Our ability to accurately retrieve the underlying inclination angle for the stars in our artificial sample allows us to place a constraint on the signal-to-noise ratio necessary to reliably extract the inclination angle. This can be placed at $\mathrm{H}/\mathrm{B}\gtrsim 4$, which can be seen in Fig.~\ref{fig: red_split_comp} as the cut where the density of points drops off greatly.

The model used for the analysis is a useful extension of the Fisher distribution, however there are a few issues that need to be addressed. The first is that there is no currently derivable analytical mean for the distribution. This can cause some issues with the interpretation of the location parameter $\mu$ of the model. In the case where $\kappa \rightarrow \infty$, $\mu$ will tend to the mean of the modified Fisher distribution, $\langle i\rangle$. For low $\kappa$ values this is not always the case due to $\mu$ not representing the mean of the distribution. We have assumed during the first part of the analysis that $\mu$ and $\langle i \rangle$ are representative of each other for large $\kappa$. This is best shown in Fig.~\ref{fig: model_diff} where the difference between $\mu$ and the mean of the distribution is shown for a range of $\kappa$ and $\mu$ values. The ranges of $\kappa$ spanned by each dataset are $20\lesssim \kappa \lesssim 145$, in the region where $\mu$ is representative of $\langle i\rangle$. This additional effect causes the hierarchical model to tend to underestimate stellar inclination angles towards zero and overestimate towards 90\mbox{$^{\circ}$}\xspace. Close to 0\mbox{$^{\circ}$}\xspace the location parameter is likely to underestimate the distribution mean whereas close to 90\mbox{$^{\circ}$}\xspace the opposite is the case (but to a much lesser degree). In an ideal case, the mean of the distribution would be the target parameter and Fig.~\ref{fig: model_diff} shows that for the majority of parameter space $\mu$ is indeed a valid alternative.

In addition to the hierarchical inference there are also assumptions involved in the fitting of the modes which should be addressed. An important assumption that goes into the asteroseismic determination of the inclination angle is that there is equipartition of energy between modes of the same azimuthal order. If this assumption were not valid then for this random sample we would expect to not see isotropy in the inferred distribution. We cannot say that the assumption has been validated, although we can certainly say that there is no evidence to suggest that the assumption is not valid.

An observational bias that needs to be considered in the determination of the inclination angle is a 90$^{\circ}$ attractor, the fact that visually it is much easier to identify a mode in a singlet (close to 0 degrees) or doublet (close to 90 degrees) configuration for a given signal-to-noise ratio than a mode in a triplet configuration. This is due to the mode visibilities and how the power is distributed amongst the $m$ components of the mode. For a given underlying height, both the singlet and doublet configuration will have larger relative heights than the triplet observed at intermediate angles. In addition, assuming an underlying isotropic distribution, observing a star with an angle close to 0 degrees is rare which results in us being much more likely to observe a star near 90 degrees. This is the main reason for selecting stars with a high \mbox{$\nu_{\mathrm{max}}$}\xspace value where mixed modes are much simpler to identify and so this effect is minimised in our sample. This will play an important role when analysing stars with much lower \mbox{$\nu_{\mathrm{max}}$}\xspace values when the mixed mode identification becomes more difficult. The location of this bias at high inclination angles is in the same location as the excess we see in both datasets close to 90\mbox{$^{\circ}$}\xspace in Fig~\ref{fig:cdfs}.

The recent work of \cite{2018arXiv180507044K} highlighted a few potential sources of bias and warned against fitting oscillation modes individually, as we have done here, since this can lead to biases in the inferred inclination angle. This issue can be seen in the Monte-Carlo estimates of our artificial datasets, whereby there are significant biases. However, we have shown that through the use of hierarchical inference it is possible to extract the inclination angle reliably even when the oscillation modes are fitted individually. Fitting all of the oscillations at the same time using a heavily parametrised model could also potentially lead to biases that a set of individual fits would not be susceptible to. An example of this is Kepler-408 where \cite{2019AJ....157..137K} resolved conflicting estimates of the inclination angle obtained from global fits to the oscillation modes with a careful treatment of the granulation background. The advantage of our method here is that our estimate of the inclination angle is unaffected by the treatment of the granulation signal due to the fact that we fit modes individually where the background is locally flat. There are alternatives such as \cite{2018MNRAS.476.1470G} who fit all the modes at once making use of a maximum likelihood estimation and a novel peak-detection scheme. Either way, the important point is that care must be taken when extracting the inclination angle of any star.

\cite{2018arXiv180507044K} also suggested that the p-dominated mixed modes should show severe blending because the ratio of the rotational splitting to the linewidth is less than 0.5. Fig.~\ref{fig: red_split_comp} shows the height-to-background ratio (HBR) as a function of the reduced splitting (which here is given by $\Gamma/\nu_{\mathrm{s}}$ rather than the reciprocal defined in \citealt{2018arXiv180507044K}) which shows that all but one of our modes across all of the data has $\Gamma/\nu_{\mathrm{s}} > 2$. This mode occurred in the artificial data set rather than the real data, and in the case of the real stars all modes lie below this limit showing that the severe blending is not an issue for these high $\nu_{\mathrm{max}}$ red-giant-branch stars.

\cite{2018A&A...616A..24G} and \cite{2018arXiv180708301M} have developed updated formalisms for describing the rotational splitting in red giant stars. The formalism presented in \cite{2018arXiv180708301M} is based upon integrating the function $\zeta$ over the frequency range between the $m=0$ and $m=\pm1$ components. The advantage of this prescription is that it can explain asymmetric rotational splittings without any need for extensions of the underlying theory. As a result, \cite{2018arXiv180708301M} reported asymmetric rotational splittings in red giant stars that exhibited fast rotation. In the case of slow rotators (such as our sample) this asymmetry is not observed. No asymmetric splittings were observed in our sample, however any small deviation from symmetry could manifest in small biases in the derived inclination angle. The degree of asymmetry however depends upon the mixing function $\zeta$ and so will be most extreme for the most g-dominated mixed modes, which in the high \mbox{$\nu_{\mathrm{max}}$}\xspace regime we are working in are rare. Therefore the possibility of observing asymmetric rotational splittings in our sample is highly unlikely. When applying this method to stars with lower \mbox{$\nu_{\mathrm{max}}$}\xspace values this will have to be an important consideration in the peakbagging and initial extraction of the inclination angle posteriors for each mode, rather than in the hierarchical method.

The possibility of asymmetric splittings is one such bias that could be present in our dataset due to assumptions made in the fitting process, such as assuming that the mode amplitudes are related by the deterministic expression given in Eqn~\ref{eqn: inc_leg}. If there are asymmetric splittings present (even slight) then the linewidths and heights of the individual azimuthal components could be affected in a similar way. Eqn~\ref{eqn: inc_leg} would therefore be correct for the underlying heights, but would then be affected by the mixing function $\zeta$. It is these sources of possible systematics that we account for by assuming an underlying inclination angle distribution for each star rather than a point estimate.

It has been suggested that searching for differences in the inclination angles inferred from p- and g-dominated mixed modes could yield evidence of core-envelope misalignment. For example, in \cite{2013Sci...342..331H} this was attempted but no such evidence was found. Unfortunately, due to the high \mbox{$\nu_{\mathrm{max}}$}\xspace values of the stars in our sample we do not have enough g-dominated mixed modes to make this comparison.

\begin{figure}
 \includegraphics[width=0.45\textwidth]{Figure9.pdf}
\caption{A scatter plot showing the distributions of the real and artificial data (drawn from the isotropic input distribution) as a function of the inverse reduced splitting ($\Gamma/\nu_{\mathrm{s}}$) and signal-to-noise (given as the ratio of the mode height to the background). The marginalised densities for each parameter are shown in the histograms coloured according to the dataset.}
\label{fig: red_split_comp}
\end{figure}

\begin{figure}
 \includegraphics[width=0.45\textwidth]{Figure10.pdf}
\caption{A contour plot showing the difference between the location parameter in the updated Fisher distribution model and the mean of the distribution.}
\label{fig: model_diff}
\end{figure}

Before asteroseismology can be used to infer full population distributions the complex selection function for the determination of the inclination angle needs to be determined. Whilst we can robustly determine the inclination angle for individual stars, there will be biases present that only become apparent when combining data at a population level, such as the relative ability to derive a given inclination angle (since some are easy to detect than others). This selection function will be dependent not only on the inclination angle, but also on the signal-to-noise ratio of the star. With the addition of a selection function then the hierarchical method presented in Section~\ref{sec: method1} can be readily extended to infer the population distribution.

\section[]{Conclusions}\label{sec:conclusions}

In this work, we have shown how asteroseismic estimates of the inclination angle of individual oscillation modes can be combined using a hierarchical Bayesian method to infer an underlying inclination angle for a single star. This has been applied to an artificial dataset that showed excellent agreement with the input inclination angles, and then applied to a sample of 123 high-$\nu_{\mathrm{max}}$ \emph{Kepler} red giants for which we report the inclination angles. The application of this method is not limited just to red-giant stars and can be readily applied to main-sequence stars to infer their inclination angles.

The hierarchical method proposed in this work could be extended to inferring the population distribution if we can derive and incorporate a selection function. This would alleviate population level biases that are not present in the individual estimates.

\section*{Acknowledgments}
We would like to thank the referee for their helpful comments and suggestions. We would also like to thank Nathalie Theme\ss l and Andr\'{e}s Garc\'{i}a Saravia Ortiz de Montellano for their very useful comments and discussions. An earlier version of this work appeared in J.S.K's PhD thesis (\url{http://etheses.bham.ac.uk/id/eprint/7658}). The research leading to the presented results has received funding from the European Research Council under the European Community's Seventh Framework Programme (FP7/2007-2013) / ERC grant agreement no. 338251 (StellarAges). This research has made use of the Exoplanet Follow-up Observation Program website, which is operated by the California Institute of Technology, under contract with the National Aeronautics and Space Administration under the Exoplanet Exploration Program. The authors acknowledge the support of the UK Science and Technology Facilities Council (STFC). Funding for the Stellar Astrophysics Centre is provided by the Danish National Research Foundation (Grant DNRF106). T. L. Campante acknowledges support from the European Union's Horizon 2020 research and innovation programme under the Marie Sk\l{}odowska-Curie grant agreement No. 792848 and from grant CIAAUP-12/2018-BPD.

This work made use of the following software: \texttt{corner} \citep{corner}, \texttt{emcee} \citep{emcee}, \texttt{Matplotlib} \citep{Hunter:2007}, \texttt{Numpy} \citep{doi:10.1109/MCSE.2011.37}, \texttt{Pandas} \citep{mckinney-proc-scipy-2010}, \texttt{scikit-learn} \citep{scikit-learn}, \texttt{Scipy} \citep{scipy}.

\bibliographystyle{mnras}


\section{introduction}

Factor models are a special kind of latent-variable models that are widely used in economics and other disciplines of research.
Due to the ubiquitous dependence across high-dimensional economic variables, it is appealing for explaining the variation of the high-dimensional economic measurements from only a small number of latent common factors.
Well-known examples include the arbitrage pricing theory \citep{Ross1976}, rank of demand systems \citep{Lewbel1991}, multiple-factor models \citep{FF1993}, components analysis in economic activities  \citep{GH1999}, and diffusion index \citep{SW2002}.
Successful applications rely heavily on the ability to correctly specify the number of factors, which is usually unknown in practice. Therefore, determining the number of factors is an essential step in applying factor models.

 One approach to this end is to utilize the eigen-structure of data matrix; see for example  \cite{F1977}, \cite{Schott1994},  \cite{YW2003}, \cite{Onatski2010} and \cite{LL2016}, amongst  others.
However, these methods require strong assumptions on the separation of eigenvalues,  and thus most of them are only applicable to fixed dimensional data.
More popular approaches are cross-validation (CV) and information criteria (IC).
{For factor model, several cross-validation methods have been proposed}; see, for example, \cite{Wold1978} and \cite{EK1982}, amongst  others.
\cite{BKSK2008} comprehensively reviewed these CV-based approaches, and concluded that most were not statistically consistent.
In addition, they recommended an alternative method based on the expectation maximization (EM) algorithm, which turned out to be extremely  computationally expensive.
As an alternative, the ICs are  more computationally effective and consistency can be guaranteed. As far as we know, \cite{CD1997} was the first paper that used ICs to determine the  number of factors.
Subsequently, other methods based on the Akaike information criterion (AIC) and Bayesian information criterion (BIC) \citep{SW1998,FLHR2000, BaiNg2002, LLS2017}  were proposed for different settings of sample size $n$ and dimension of variable, $p$.
When the variance of the idiosyncratic component is relatively small, these methods are able to provide very efficient estimates of the number of factors \citep{BaiNg2002}.
However, the IC-based approaches are usually not fully data-driven since their penalties depend on predetermined tuning parameters.
Even though several guides are proposed to mitigate the effect of  tuning parameters, the performances of these methods are not stable when the signal-to-noise ratio is relatively small \citep{Onatski2010}, i.e., the variation of the idiosyncratic component is relatively large, or when the idiosyncratic component has heavy tails in distribution.
Unfortunately, both of which are common in financial data.

In this paper, we propose to estimate the number of factors by a computationally-efficient CV method.
Note that the conventional CV methods, e.g., those used in the linear regression, only leave one or several \emph{observations} out.
In contrast, factor models involve a matrix, where both observations (or rows) and variables (or columns) play similar roles mathematically in the analysis.

Thus, we propose a double cross-validation (DCV) method that applies CV first to the \emph{observations} and then to the \emph{variables} in the {observations}.
Theoretically, we show that the method is consistent for high-dimensional data and allows dependency between the idiosyncratic errors.
The method  is thus applicable to approximate factor models.
In addition, computationally, the second CV can be easily calculated by the linear regression error using the full data \citep{Shao1993}, while the first CV can be done in a $K$-fold manner.
Simulation studies show that the proposed approach performs satisfactorily even when the idiosyncratic  error has relatively big variance or homoscedasticity, or heavy tail,  which are especially relevant for economic data.

The rest of this paper is organized as follows. Section~\ref{sec:method} reviews the basic notation of factor analysis.
Section~\ref{DCVsection} presents the DCV approach.
Section~\ref{sec:theory} establishes statistical consistency of the approach, while the proofs of those properties are given in the Appendix.
Sections~\ref{sec:simulation}--\ref{sec:application} present a set of simulation studies and  an empirical application to assess the finite sample performance of DCV. Concluding remarks are given in Section~\ref{sec:discussion}.

\section{Factor model and its  factors and  loadings}\label{sec:method}

For a $p$-dimensional random vector ${\bm x}=(X_1, X_2,...,X_p)^\top $, the factor model assumes the following \emph{bilinear} structure
\begin{equation*}
{\bm x}={\bm L^0{\bm f}^0}+{\bm e},
\end{equation*}
where ${\bm L}^0\in\mathbb{R}^{p\times d_0}$ is the loading matrix, ${\bm f}^0=(f_{1}, f_{2}, ... , f_{d_0})^\top $ is  the vector of $ \ d_0 \  $ common factors. Together ${\bm L^0{\bm f}^0}$ is called the common component of the data,  and ${\bm e}=(e_{1},e_{2},...,e_{p})^\top $ the idiosyncratic component. Let $\boldsymbol\Sigma_{\bm e}=\mathbb{E}[\bm e\bm e^\top]$, $\boldsymbol\Sigma_{\bm f_0}=\mathbb{E}[{\bm f_0\bm f_0^\top}]$, and $\boldsymbol\Sigma_{\bm x}=\mathbb{E}[\bm x\bm x^\top]$. Under this model, the variation of the $p$-dimensional variable is mainly generated by a small number of common factors.
In the statistical literature, the idiosyncratic component are usually assumed to be independent of the common component, and have diagonal covariance matrix, in which case {$\boldsymbol\Sigma_{\bm x}={\bm L}^0\boldsymbol\Sigma_{\bm f_0}{{\bm L}^{0}}^\top +\boldsymbol\Sigma_{\bm e}$}. However, such independence assumption and diagonal structure are usually not realistic in financial data or macroeconomic data where factor analysis are often applied. In this paper, we consider the approximate factor model where the idiosyncratic components can be weakly dependent.

Let  ${\bm X}=({\bm x}^\top _1,{\bm x}^\top _2,...,{\bm x}^\top _n)^\top $ be random samples of  $n$ observations on $ p$ variables. We have the following matrix form for the factor model:
\begin{equation}  \label{FM}
{\bm X}={\bm F}^0{{\bm L}^{0}}^\top +{\bm E},
\end{equation}
with ${\bm F}^0=({\bm f}_1^0,{\bm f}_2^0,...,{\bm f}_n^0)^\top $ and  ${\bm E}=({\bm e}^\top _1,{\bm e}^\top _2,...,{\bm e}^\top _n)^\top $.
According to the eigen-decomposition, ${\bm X^\top \bm X}$ has the  expression
\begin{equation*}
{\bm X^\top \bm X}=\boldsymbol\Phi\boldsymbol\Lambda\boldsymbol\Phi^\top,
\end{equation*}
where $\boldsymbol\Lambda$ is a $p\times p$  diagonal matrix  with diagonal  entries being  the eigenvalues of $\bm X^\top\bm X$, $\boldsymbol\Phi=(\boldsymbol\phi_1,\boldsymbol\phi_2,...,\boldsymbol\phi_n)\in\mathbb{R}^{n\times n}$ is the eigenvector matrix of ${\bm X^\top \bm X}$ such that $\boldsymbol\Phi^\top\boldsymbol\Phi={\bm I}_p$. Then, with working number of factors $d$, the estimators of factor loadings and factors  are respectively
\begin{equation}\label{FL1}
\widetilde{\bm L}^d=\sqrt{p}(\boldsymbol\phi_1,\boldsymbol\phi_2,...,\boldsymbol\phi_d) \ \  \mbox{and} \ \  \widetilde{\bm F}^d=\frac{1}{{p}}{\bm X}\widetilde{\bm L}^d.
\end{equation}
Denote by $\text{Tr}(\bm A)$ the trace of matrix $\bm A$. It is easy to see that
\begin{equation*} \label{LF}
(\widetilde{\bm F}^{d},\widetilde{\bm L}^{d})=\underset{\substack{{\bm F}^d\in\mathbb{R}^{n\times d},{\bm L}^d\in\mathbb{R}^{p\times d}\\{{\bm L}^{d}}^\top {\bm L}^d/p={\bm I}_d}}{\text{argmin}}\{\text{Tr}[({\bm X}-{\bm F}^d{{\bm L}^{d}}^\top )({\bm X}-{\bm F}^d{{\bm L}^{d}}^\top )^\top ]\}.
\end{equation*}

\section{Double Cross Validation for Number of Factors}\label{DCVsection}

Our basic idea to specify the number of factor $d_0$, is based on the prediction error in a cross-validated manner. Write the model element-wise as
\begin{equation}\label{reg1}
{x}_{is}={{\bm f}_i^0}^\top {\bm l}^0_s+{{e}}_{is}, \ \ \ \ i=1, ..., n, \ \  s=1,2,...,p.
\end{equation}
CV needs to make prediction of $ {x}_{is} $ based on other elements except for itself.

Because neither $ {{\bm f}_i^0} $ nor $  {\bm l}^0_s $ in \eqref{reg1} is observable, we implement the prediction by two-stage fitting: leave-\emph{observation}-out  and leave-\emph{variable}-out.

The first stage is to estimate $  {\bm l}^0_s $ from the data by the $K$-fold CV. Divide the rows of $ {\bm X} $ into $K$ folds, $ 1< K \le n  $. Denote them by $M_1, ..., M_K$. Let $n_k=\#M_k$ be the number of elements in $M_k$, and ${\bm X}_{-M_k}$ be the sub-matrix of $ {\bm X} $ with  rows in $ M_k $ being removed. Apply  the eigen-decomposition to ${\bm X}_{-M_k}$ to obtain corresponding matrices $\widetilde{\bm F}^{k,d}$ and $\widetilde{\bm L}^{k,d}$ as in \eqref{FL1}.
Here, we use the superscript to highlight the fact that the estimator is obtained from the data with the rows in $M_k$ removed, with working number of factors $d$. Similar to \cite{BaiNg2002}, we rescale the estimated factor loading matrix and let
\begin{equation*}
\widehat{\bm L}^{k,d}= (\frac{1}{n-n_k}{\bm X}^\top_{-M_k}{{\bm X}_{-M_k}}) \widetilde{\bm L}^{k,d} = (\widehat{\bm l}^{k,d}_1, ..., \widehat{\bm l}^{k,d}_p).
\end{equation*}

In the second stage, we replace $ {\bm l}^0_s $ in \eqref{reg1} by $ \widehat{\bm l}^{k,d}_s $ and rewrite \eqref{reg1} as a regression model
\begin{equation}\label{reg2}
{x}_{is}={{\bm f}_i^0}^\top \widehat{\bm l}^{k,d}_s + {{e}}'_{is}, \ \ \ \ i \in M_k, \ \  s=1,2,...,p,
\end{equation}
where $ e'_{is} $ is the regression error in place of $ e_{is} $ due to the replacement. This time,  $\widehat{\bm l}^{k,d}_s, s=1,2,...,p$ are known and treated as the `regressors', but ${\bm f}^0_i$ is treated as  `regression coefficients' that need to be estimated.
By leaving $x_{is}$ out,  we estimate ${\bm f}^0_i$ by
\begin{equation*}
 \widehat {{\bm f}}_{i,s}^{k,d} = \underset{\bm f}{\text{argmin}} \sum_{t\neq s, t=1}^p ( {x}_{it} -\bm f^\top \widehat{\bm l}^{k,d}_t)^2 ,\  i \in M_k.
\end{equation*}
Thus we can predict $x_{is}$ by $(\widehat {{\bm f}}_{i,s}^{k,d})^\top \widehat{\bm l}^{k,d}_s $.
The average squared prediction error for $\bm x_i = (x_{i1}, ..., x_{ip})^\top$ is
\begin{equation*}
 V_i^{k,d} = \frac{1}{p}\sum\limits_{s=1}^p [x_{is} - (\widehat {{\bm f}}_{i,s}^{k,d})^\top \widehat{\bm l}^{k,d}_s]^2, \ \  i \in M_k.
\end{equation*}
The calculation of $  V_i^{k,d} $ can be much simplified  as shown in \cite{Shao1993}, i.e.,
\begin{equation*}
V_i^{k,d} =\frac{1}{p}\sum\limits_{s=1}^p(1-w^{k,d}_s)^{-2}(x_{is}-(\widehat{\bm f}_{i}^{k,d})^\top\widehat{\bm l}^{k,d}_s)^2,
\end{equation*}
where  $\widehat{\bm f}_i^{k,d}=(\widehat{\bm L}^{k,d\top}\widehat{\bm L}^{k,d})^{-1}\widehat{\bm L}^{k,d\top}\bm x_i$ is the conventional least squares estimator for \eqref{reg2} and $w^{k,d}_s$  is the $s$-th diagonal element of projection matrix $\mathcal{P}_{\widehat{\bm L}^{k,d}}=\widehat{\bm L}^{k,d}(\widehat{\bm L}^{k,d\top}\widehat{\bm L}^{k,d})^{-1} \times \widehat{\bm L}^{k,d\top}$.
Finally, consider the averaged  prediction error over all the elements
\begin{equation*}
DCV(d)=\frac{1}{n} \sum\limits_{k= 1}^K \sum\limits_{i \in M_k}^nV_i^{k,d}.
\end{equation*}
 Let $  d_{\max}$ be a fixed positive integer that is large enough such that  $  p> d_{\max} > d_0$.  The DCV estimator for the number of factors is given by
\begin{equation}\label{definek}
\widehat{d}_{DCV}=\underset{0\leq d\leq d_{\max}}{\text{argmin}}DCV(d).
\end{equation}

The above approach has similarity with the vanilla row-wise CV \citep{BKSK2008}, which also predicts ${\bm x}_{i}$  by two steps. Its first step is exactly the same as ours. However, in the second step, row-wise CV uses the full data instead of cross validation method,  leading to overfitting and inconsistency of the method. To fix this problem, \cite{EK1982}  suggested the element-wise cross validation. Although their methods solve the problem of overfitting, they are costly and involve immense computation.

Our DCV  offers a simpler solution. By utilizing CV in both the first stage of leaving ``observations'' out, and the second regression step of leaving   ``variables'' out, our method  guarantees that the prediction of each element does not use any information from itself, and thus the consistency is ensured as shown in the next section. In addition, the implement of $K$-fold CV can significantly reduce the computational complexity, especially when the number of rows is large. It is interesting to see that the $K$-fold CV in the first stage will also help in simplifying the calculation in the second stage. Because of this appealing nature, to further facilitate the calculation, we can transpose $ {\bm X} $ when $n<p$.

\section{Asymptotic consistency of the estimation}\label{sec:theory}
In this section, we investigate the consistency of our derived estimator.  We start with the following assumptions.
\begin{ass}\label{ass1}
The eigenvalues of ${\bm F}^{0\top}{\bm F}^0/n$ are bounded away from zero and infinity, and  that $\mathbb{E}\|{\bm f}^0_i\|^4<\infty$ for $ i=1, ..., n$.
\end{ass}
\begin{ass}\label{ass2} The eigenvalues of ${\bm L}^{0\top}{\bm L}^0/p$ are bounded away from zero and infinity, and that $\mathbb{E}\|{\bm l}^0_i\|^4<\infty$ for  $i=1, ..., p$.
 \end{ass}

\begin{ass}\label{ass3} There exists a constant $M<\infty$ such that
\begin{enumerate}
\item[1.]  $\mathbb{E} e_{it}=0$  and  $\mathbb{E} e_{it}^4\leq M$ for $1\leq i\leq n$ and $1\leq t\leq p$;

  \item[2.] $\gamma_p(i,j) = \mathbb{E}({\bm e}_i^\top{\bm e}_j/p)\leq M$  and $\frac{1}n\sum\limits_{j=1}^n\gamma^2_p(i,j)\leq M$ for all $i$;

 \item[3.] $\mathbb{E}({\bm e}_{is}{\bm e}_{it})=\tau_{st,i}$ with $\tau_{st,i}\leq |\tau_{st}|$ for some $\tau_{st}$ and for all $i$; in addition,  $\frac{1}p\sum\limits_{s=1}^p\sum\limits_{t=1}^p\tau_{st}$ $\leq M$;

 \item[4.] $\mathbb{E}({\bm e}_{is}{\bm e}_{jt})=\tau_{ij,st}$ and $\frac{1}{np}\sum\limits_{i=1}^n\sum\limits_{j=1}^n\sum\limits_{s=1}^p\sum\limits_{t=1}^p\tau_{ij,st}\leq M$;

\item[5.]  $\frac{1}{  n }\mathbb{E}|\sum\limits_{j=1}^n[e_{jt}e_{js}-\mathbb{E}(e_{jt}e_{js})]|^2\leq M$ for all $(s,t)$;

\item[6.]  $\|\frac{1}{\sqrt{n}}\sum\limits_{i=1}^n{\bm f}^{0}_{i}e_{is}\|^2\leq M$ for $1\leq s\leq p$;

\item[7.]  $\|\frac{1}{\sqrt{p}}\sum\limits_{s=1}^p{\bm l}^{0}_{s}e_{is}\|^2\leq M$ for $1\leq i\leq n$.

\end{enumerate}
\end{ass}

Assumptions \ref{ass1}--\ref{ass3} are commonly used in the  analysis of  factor models, for example, in \citet{BaiNg2002} and \citet{LLS2017}.  Assumptions \ref{ass1} and  \ref{ass2} together ensure that each factor plays a nontrivial role in contributing to the variation of ${\bm X}$. Unlike the strict factor model that assumes all entries of $\bm E$ are I.I.D., Assumption~\ref{ass3} allows weak dependency for elements of  $\bm E$, making our methods applicable to the approximate factor models.  These assumptions also indicate that  all the  eigenvalues of the covariance matrix of common components would dominate the eigenvalues corresponding to idiosyncratic components, which is crucial to the identifiability for the approximate factor models.

\begin{ass}\label{ass4}
$\sup\limits_k n_k=o(n^{1/3})$ and
 $K\sup\limits_kn_k/n< \infty$, where $n_k=\#M_k$ is the number of elements in $M_k$.

\end{ass}
Assumption \ref{ass4} is a weak condition on $K$-fold CV in the first stage. The number of elements in each fold can either be fixed or tend to infinity. In particular, this assumption includes the leave-one-out CV.
\edef\oldass{\the\numexpr\value{ass}+1}

\setcounter{ass}{0}
\renewcommand{\theass}{\oldass.\alph{ass}}

\begin{ass}\label{ass5a}  The idiosyncratic components satisfies
\begin{enumerate}
\item[1.] $\mathbb{E}[e_{it}^2\mid {\bm F}^0_{-M_k}{\bm L}^{0\top}]=\sigma^2$ for $i\in M_k$, where ${\bm F}^0_{-M_k}$ is a submatrix of ${\bm F}^0 $ with rows in $M_k$ being removed;

\item[2.] $\lambda_1({\bm E^\top}{\bm E}/n)<2 \sigma^2$.
\end{enumerate}

\end{ass}

\begin{ass}\label{ass5b} Let $p/n\to \rho\in(0,\infty)$.  There exists a $p\times p$ positive definite matrix $\boldsymbol\Sigma_p$ such that ${\bm e}_i=\boldsymbol\Sigma_p^{1/2}{\bm y}_i$, where ${\bm y}_i=( y_{i1},y_{i2},...,y_{ip})^\top$,
\begin{enumerate}
\item[1.] $\mathbb{E} y_{it}=0$, $\mathbb{E} y_{it}^4<\infty$,  and $y_{it}$ are I.I.D. for $1\leq i\leq n, 1\leq t\leq p$;

\item[2.]  ${\bm y}_i$ is independent of ${\bm F}^0_{-M_k}$;

 \item[3.] The spectral distribution of $\boldsymbol\Sigma_p$ is convergent, i.e., there exists a distribution function $H$ such that $\frac{1}{p}\sum\limits_{s=1}^p{\mathds{1}}(\lambda_s(\boldsymbol\Sigma_p)\leq t)\to H(t)$ as $p\to \infty$;

\item[4.] $2\text{\normalfont{diag}}(\boldsymbol\Sigma_p) - \boldsymbol\Sigma_p $ is  a positive definite matrix, where $\text{\normalfont{diag}}(\boldsymbol\Sigma_p)$ is a diagonal matrix that has the same diagonal elements as $\boldsymbol\Sigma_p$.
\end{enumerate}
\end{ass}
\let\theass\origtheass

Assumptions \ref{ass5a} and  \ref{ass5b} are two technical assumptions, either of which ensures the viability of the proposed method.  They also allow weak dependence among idiosyncratic components as well as the observations, which generalize the I.I.D. assumption required in the strict factor model.  In comparison, \cite{BaiNg2002} needs weaker assumptions on  the dependence than ours, possibly due to the relatively large penalty in their methods. Of course, this large penalty is at the cost of selecting inappropriately fewer factors as shown in our simulation study.

\begin{thm} \label{thm1}
Suppose  Assumptions \ref{ass1}--\ref{ass4} hold. If  either Assumption \ref{ass5a}  or \ref{ass5b} holds, then the DCV estimator \eqref{definek} is consistent as both $ n $ and $ p $ tend to infinity, i.e.,
\begin{equation*}
\lim\limits_{n,p\to\infty}Prob(\widehat{d}_{DCV}=d_0)=1.
\end{equation*}
\end{thm}

{\color{black}
\begin{rem}  Recent empirical findings suggest that the factor structure of economic data may change with  sample size and dimensionality of variable \citep{Onatski2010,  JLN2015,LLS2017}. In fact, Theorem \ref{thm1} can be easily extended to the case that $ d_0 $ varies slowly with both $ n $ and $ p $. {\color{red}DCV} is thus applicable to practical data with changing structure of factorization.
\end{rem}}

\section{Simulation Study}\label{sec:simulation}
In this section, we conduct simulation studies to compare the finite sample performances of the DCV with the other methods, including
the {panel} criterion $IC_1$ \citep{BaiNg2002} and $Ladle$ \citep{LL2016}.
As shown in \cite{LL2016},  $Ladle$ estimator outperforms all the other eigen-structrue based methods in determining the number of factors, and is thus used  as a representative of those methods.
Similarly, for the IC-based methods, \cite{BaiNg2002} showed that their criterion outperforms other information criteria such as $AIC$ and $BIC$.
We do not consider other cross-validation methods as they are either inconsistent or suffer from excessive computational burden.
For our $DCV$, we consider both leave-one-out CV ($DCV_1$) and 10-fold CV ($DCV_{10}$) in the first stage.
Matlab codes for the calculations are available at \href{https://github.com/XianliZeng/Double-Cross-Validation}{https://github.com/XianliZeng/Double-Cross-Validation}.

\begin{figure*}[ht!]
\vspace{-0.3cm}
\centerline{\includegraphics[width=1.0\linewidth]{model1.eps} } \vspace{-0.5cm}
 \caption{\footnotesize The relative frequencies of correctly selecting the number of factors for model~\eqref{SM1} with independent Gaussian idiosyncratic errors (E1). In each panel, we use red line for $DCV$, black line for $DCV_{10}$, {blue} line for $IC_1$ and  green line for \textit{Ladle}.\vspace{-0.5cm}
}
\label{fig1}
\end{figure*}
\begin{figure*}[ht!]
\vspace{-0.3cm}
\centerline{\includegraphics[width=1.0\linewidth]{model2.eps} } \vspace{-0.5cm}
 \caption{\footnotesize The relative frequencies of correctly selecting the number of factors for  model \eqref{SM1} with independent $t$-distributed idiosyncratic errors (E2). In each panel, we use red line for $DCV$, black line for $DCV_{10}$, blue line for $IC_1$ and  green line for \textit{Ladle}.\vspace{-0.5cm}
}
\label{fig2}
\end{figure*}
\begin{figure*}[ht!]
\vspace{-0.3cm}
\centerline{\includegraphics[width=1.0\linewidth]{model3.eps} } \vspace{-0.5cm}
 \caption{\footnotesize The relative frequencies of correctly selecting the number of factors for   model~\eqref{SM1} with heteroskedastic idiosyncratic errors (E3). In each panel, we use red line for $DCV$, black line for $DCV_{10}$, blue line for $IC_1$ and  green line for \textit{Ladle}.\vspace{-0.5cm}
}
\label{fig3}
\end{figure*}

\begin{figure*}[ht!]
\vspace{-0.3cm}
\centerline{\includegraphics[width=1.0\linewidth]{model4.eps} } \vspace{-0.5cm}
 \caption{\footnotesize The relative frequencies of correctly selecting the number of factors for  model~\eqref{SM1} with serial correlated  idiosyncratic errors (E4). In each panel, we use red line for $DCV$, black line for $DCV_{10}$, blue line for $IC_1$ and  green line for \textit{Ladle}.\vspace{-0.5cm}
}
\label{fig4}
\end{figure*}

\begin{figure*}[ht!]
\vspace{-0.3cm}
\centerline{\includegraphics[width=1.0\linewidth]{model5.eps} } \vspace{-0.5cm}
 \caption{\footnotesize The relative frequencies of correctly selecting the number of factors for  model~\eqref{SM1} with cross-sectional correlated idiosyncratic errors (E5). In each panel, we use red line for $DCV$, black line for $DCV_{10}$, blue line for $IC_1$ and  green line for \textit{Ladle}.\vspace{-0.5cm}
}
\label{fig5}
\end{figure*}

We simulate data from a factor model with $d_0 = 5$:
\begin{equation}
 x_{is}=\sum_{j=1}^{5}f_{ij}l_{sj}+\sqrt{\theta}e_{is}; \quad f_{ij},l_{sj}\sim N(0,1)  \ \ i=1, ..., n, \ \  s=1,2,...,p.  \label{SM1}
\end{equation}
For the idiosyncratic errors $ e_{is} $, the following settings are considered:
\begin{description}
\item[] (E1) independent Gaussian:  $ e_{is}\sim N(0,1)$;
\item[] (E2) independent t-distributed:  $ e_{is}\sim t_3 $;
\item[] (E3) heteroskedastic:  $ e_{is} \sim N(0,\delta_s), $ $ \delta_s = 1 $ if $s$ is odd,  and 2 if $s$ is even;
\item[] (E4) serial correlated: $ e_{is}=0.3 e_{i,s-1}+\nu_{is} $,  $\nu_{is}\sim N(0,1)$;
\item[] (E5) cross-sectional correlated: $  e_{is}=\sum\limits_{j=-10}^{10} 0.15^j \nu_{i-j,s}$,  $\nu_{is}\sim N(0,1)$.
\end{description}

In all simulations, we select $d$ from $\{1,2,...,d_{\max}\}$ where $d_{\max}= 8 $ as in \cite{BaiNg2002}.
We consider the combinations of { $n$ and $p$} with $n=\{40,160,640\}$ and $p=\{30,90, 270\}$.
To see the effect of  the signal-to-noise ratio on the methods, $\theta$ varies in a range of intervals depending on different settings of idiosyncratic errors and different $ n $ and $p$.
These ranges are made clear in each panel of the figures below.
All comparisons  were evaluated based on 1000 replications.
Figures~\ref{fig1} to~\ref{fig5} summarize the simulation results, where the curves represent the frequencies of selecting the right number of factors at different noise level $\theta$.

Generally, the simulation results suggest that the performance of all the methods becomes better when $n$ and $p$ get larger, which lends support to the asymptotic consistency of the methods.
$Ladle$ usually shows better performance than $IC_1$, possibly due to its utilization of  both eigenvalues and eigenvectors.
Notably, $DCV_1$ and $DCV_{10}$ stand out and perform the best in most scenarios.
They possess a much higher frequency of correct selection when $\theta$ is large, implying their ability to select the right number of factors even when the {factors' contribution to the variation of the variables,  i.e., signal-to-noise ratio,} is low, which is particularly the situation in economic data.
Moreover, $DCV_1$ and $DCV_{10}$ present almost
identical performance for all the models, implying that the performance of the DCV is not sensitive to the choice of~$K$.

For different settings of the idiosyncratic errors, we have some additional observations.
Figure \ref{fig2} shows that $IC_1 $ is not so robust to the data with heavy-tail especially when the size of the data is small.
Figure \ref{fig3} implies that $IC_1 $ is also affected adversely by the heteroskedastic errors.
Note that both heavy tail and  heteroskedasticity are stylized facts in financial data.
In contrast, for these two types of data,  $DCV_1$ and $DCV_{10}$ still demonstrate very stable and efficient ability in selecting the number of factors.
In Figures \ref{fig4} and \ref{fig5}, where the idiosyncratic errors have serial or cross-sectional correlation,  DCV has inferior performance than $IC_1$ and $ Ladle $ when the variance of idiosyncratic components is small and  both $ n $ and $ p $ are small.
However, DCVs quickly gains advantage when either $ n $ or $ p $ increases.

\section{An empirical application: Fama-French  three-factor model}\label{sec:application}
\begin{figure*}[p]
\vspace{-0.5cm}
\centerline{Validation of three-factor model  based on value weighted return of 25 portfolios}
\centerline{\includegraphics[width=1\linewidth, height=0.56\linewidth]{figreal1.pdf}}

\centerline{Validation of three-factor model  based on value weighted return of 100 portfolios}
\centerline{\includegraphics[width=1\linewidth, height=0.56\linewidth]{figreal2.pdf}} \vspace{-0.5cm}
 \caption{\footnotesize The  frequencies of selected numbers of factors for the value weighted returns of 25 and 100 portfolios. Panels in each column represents one method; panels in the same rows are based on the same data in fixed period of years.
}
\label{fig6}
\end{figure*}

In asset pricing and portfolio management, the so-called three-factor model of \cite{FF1993} is widely applied to describe the stock returns.
It shows that the excess return of a stock or a portfolio ($R_{it}-{r_{ft}}$) can be satisfactorily explained by three common factors:  (1) excess return of the market portfolio ($R_{Mt}-r_{ft}$), (2) the outperformance of small-cap companies versus big-cap companies ($SMB$), and (3) the outperformance of companies with  high book to market ratio versus those with small book to market ratio ($HML$), i.e.,
\begin{equation*}
R_{it}-r_{ft}=\beta_1(R_{Mt}-r_{ft})+\beta_2SMB_t+\beta_{3}HML_t+\epsilon_{it},
\end{equation*}
where $ r_{ft} $ is the risk-free return.
Next,  we use the data provided by Professor  Kenneth R. French\footnote{http://mba.tuck.dartmouth.edu/pages/faculty/ken.french/data\_library.html} to mutually validate our methods and the three-factor model.
Stocks are divided into 6, 25 or 100 groups according to company's capitalization (market equity) and  Book-to-Market ratio, and portfolios are constructed with  value weighted daily returns or equally weighted daily returns in each group.
As those portfolios were compiled based on company's data at the end of June of every year, we set the first of July of every year as a starting date of a period.
It is understandable that the  portfolios change from year to year, thus we only consider periods of 1 to 3 years in our calculation.

Figure \ref{fig6} presents the frequencies of the numbers of factors detected by the $DCV_1$, $DCV_{10}$, $IC_1$ and $Ladle$ in two of the data sets with $d_{\max}$ set at $15$.
In most situations, all the  approaches identify three common factors for the returns of portfolios, which are in line with the \emph{three}-factor model of \cite{FF1993}.
In general,  DCV selects three factors most often amongst the approaches, showing its superior performance in practice.

\section{Concluding Remarks}\label{sec:discussion}
In the literature, both CV and IC based approaches were popularly used for determining the complexity of a model, such as  the number of factors.
As noticed in \cite{BKSK2008}, most CV approaches for factor models are not consistent due to their insufficient validation, while existing remedies for the consistency is very computationally expensive.
In contrast, ICs are consistent and easy to implement, but they depend on  predetermined penalty  functions \citep{BaiNg2002}.
Simulation results not reported here show that the penalty function in \cite{BaiNg2002} could  be modified to be more adaptive to the data.
In addition, ICs are unstable when the signal-to-noise ratio is relatively low \citep{Onatski2010}.

By validating the model twice, the proposed DCV method does not only ensure the consistency under mild conditions, but is also easy to implement. %
Because DCV is based on prediction error, it automatically selects the number of factors that balances the model complexity and stability.
Our simulation studies also demonstrate its superior efficiency over the existing methods most of the time.
The only exception is when there is strong serial dependence in the idiosyncratic errors.
However, this deficiency disappears when the sample size and dimension of the variables increase.
The advantages of our method are more  pronounced for data with relative large variation, heteroscedasticity  or heavy tails in the idiosyncratic errors. The method is thus particularly relevant for financial data.

\bibliographystyle{apalike}