Past exam of the mathematics course of the University of Cambridge / 2007 / iii / Paper 48 (Time Series and Monte Carlo Inference)

Define the autoregressive polynomial $\Phi (z)=1-{\sum }_{k=1}^p{\varphi }_k{z}^k$ and the moving-average polynomial $\Theta (z)=1+{\sum }_{k=1}^q{\theta }_k{z}^k$. The standard causality and invertibility root criteria for an ARMA model areIn the usual minimal autoregressive moving-average model, the two polynomials have no common roots. Under that convention, these are the necessary and sufficient root conditions for a causal and invertible representation with the given noise. Without minimality they remain sufficient; shared factors must be considered before claiming necessity.

To see why the conditions work, both ratios $C(z)=\Theta (z)/\Phi (z)$ and $D(z)=\Phi (z)/\Theta (z)$ are analytic on a disc of radius greater than one. Their power series coefficients therefore decay geometrically and are absolutely summable. The backshift operator $B{X}_t=X_{t-1}$ then gives a well-defined mean-square convergent causal time-series representationFiltering white noise by these fixed coefficients gives a weakly stationary process. The identity $\Phi (B)C(B)=\Theta (B)$ verifies the original equation. Multiplication by the absolutely summable inverse of $\Phi (B)$ also shows uniqueness among finite-variance stationary solutions. Similarly $D(B)X_t={ϵ}_t$ gives the invertible time-series representation.

Comparing coefficients in $\Phi (z)C(z)=\Theta (z)$ derives the first of the ARMA coefficient recursions:Here $c_j=0$ for negative $j$, ${\theta }_0=1$, and ${\theta }_j=0$ for $j>q$, so in particular $c_0=1$. For the inverse filter, compare coefficients in $\Theta (z)D(z)=\Phi (z)$. With $a_0=1$, $a_j=-{\varphi }_j$ for $1\le j\le p$, $a_j=0$ for $j>p$, and $d_j=0$ for negative indices, this givesThe minus sign in the inverse recursion comes from moving the nonconstant moving-average terms to the other side.

For the first-order model, the standard uncancelled conditions are $|\varphi |<1$ and $|\theta |<1$. The geometric power series expansions giveThus the ARMA(1,1) causal and inverse coefficients areAt $j=1$ these give $c_1=\varphi +\theta$ and $d_1=-(\varphi +\theta )$, including zero-valued parameters without any negative powers. If $\theta =-\varphi$, the filtered process reduces to $X_t={ϵ}_t$ and all later coefficients vanish. This common factor cancellation in an ARMA model is why the root inequalities should not be advertised as necessary for every nonminimal parameterization. The stated inequalities still ensure uniqueness for the original equation as well as the desired causal and invertible filters.

The scalar process has zero mean, because white noise has zero mean. Expanding products and retaining only matching noise indices gives its autocovariance function:The finite second moment and covariance depending only on the lag prove weak stationarity. With the two-sided angular-frequency convention, the spectral density of a stationary process isIt is nonnegative because the bracket is $|1+\theta e^{i\lambda }{|}^2$, and it is even because the cosine is even. Hence $f(-\lambda )=f(\lambda )$. No restriction on the real parameter $\theta$ is needed for stationarity here; invertibility is a separate condition.

For the bivariate calculation let $a={\sigma }_u^2$, $b={\sigma }_v^2$, $D=\mathrm{diag}(a,b)$, and $\Theta =({\theta }_{ij})$. Uncorrelated white noise processes here means zero cross-covariances at every pair of times. Therefore the vector noise hasThe filtered vector is the multivariate moving-average model $Z_t=W_t+\Theta W_{t-1}$. Linearity gives $\mathbb{E}{X}_t=\mathbb{E}{Y}_t=0$. Put $\Gamma (h)=\mathbb{E}[Z_t{Z}_{t+h}^T]$, with exactly the orientation of the printed covariance definition. At lag zero, the present and past noise terms are uncorrelated, giving $D+\Theta D{\Theta }^T$. At lag one the only shared noise is $W_t$, contributing $D{\Theta }^T$; at lag minus one it contributes $\Theta D$. Thus the covariance and cross-spectrum of a vector MA(1) are determined byExplicitly,In particular the off-diagonal positive- and negative-lag covariances need not coincide. Reading the upper-right entries and taking their Fourier series gives the requested cross-spectrum:Every coefficient is real, so replacing $\lambda$ by $-\lambda$ conjugates both complex exponentials. Therefore $f_{XY}(-\lambda )=\overline{f_{XY}(\lambda )}$. Unlike the scalar auto-spectrum, this cross-spectrum may be complex and need not be even: its imaginary part is $(b{\theta }_{12}-a{\theta }_{21})\sin \lambda /(2\pi )$.

For $0<u<1$ define the generalized inverse, or quantile function, byIt is finite because the cumulative distribution function tends to zero and one at the two infinities. Monotonicity and right continuity of $F$ give the key equivalenceFor the forward implication, choose points in the defining set approaching its infimum from above. Right continuity gives $F(Q(u))\ge u$, and then monotonicity gives $F(x)\ge u$ whenever $Q(u)\le x$. The reverse implication follows immediately because $x$ then belongs to the defining set. Consequently, for a continuous uniform distribution draw $U$,This proves inverse transform sampling for discontinuous as well as continuous distribution functions.

There is a minor endpoint defect in the printed codomain. The true quantile endpoint at probability one is $Q(1)=\inf \{x:F(x)=1\}$, which equals $+\infty$ for an unbounded-above normal distribution. Thus it need not be real at one. Use the extended-real value there, or assign an arbitrary finite value at this single probability level. Since $P(U=1)=0$, either choice gives the same sampled distribution; the latter supplies a real-valued version on all of $(0,1]$. The endpoint has no effect on the inverse-transform law.

For rejection sampling, choose a proposal probability density function $g$ that can be sampled and a finite constant $M$ with $f(y)\le Mg(y)$ wherever the target has positive mass. Generate independent pairs $(Y_j,U_j)$, with $Y_j$ of density $g$ and $U_j$ uniform. Accept $Y_j$ if $U_j\le f(Y_j)/(Mg(Y_j))$; otherwise discard the pair and repeat with fresh independent draws. At points where both densities vanish the acceptance rule may be defined arbitrarily. For any measurable set $A$,Thus one trial accepts with probability $1/M$, and the conditional density given acceptance is $f$. More explicitly, summing over the index of the first accepted proposal givesThe probability of never accepting is zero, since $M<\infty$, and the expected number of proposals is $M$. The algorithm terminates almost surely and samples the target density exactly.

Substitute $x={\sin }^{2}v$, with $0<v<\pi /2$, when integrating the target density. Since $dx=2\sin v\cos vdv$, its cumulative distribution function isInverting gives the arcsine distribution samplerThe inverse transform sampling proof in the root solution shows that it has the required probability density function. Use successive independent uniforms for successive samples. A draw at the endpoint one can be discarded if values strictly inside the support are desired; this event has probability zero under the ideal continuous uniform distribution.

Group the two terms containing $p$ to obtain the mixture distributionThe first component is a Standard Cauchy distribution. The second is normalized becauseFor this symmetric component, let $R=|X|$. Its density is $2r/(1+r^2{)}^2$, and hence its cumulative distribution function and quantile function areThe independent sign is positive or negative with equal probability. For one sample, use the next three independent uniforms $U,V,H$ in the given sequence and returnThe Cauchy branch follows from its quantile function, while the other branch follows from the magnitude law and fair sign just derived. Independence of the branch selection and component draws therefore gives exactly the displayed mixture distribution. Discard the zero-probability endpoint draw $V=1$ to keep the computation finite, and use a fresh block of three uniforms for the next sample. This is the symmetric rational mixture sampler.

One can alternatively use the root solution's rejection sampling with a Standard Cauchy distribution proposal. The density ratio is $p+\pi (1-p)|x|/(1+x^2)$, bounded by $M=p+\pi (1-p)/2$ because $2|x|\le 1+x^2$. Accept the Cauchy draw $Y$ with probability $[p+\pi (1-p)|Y|/(1+Y^2)]/M$ using a fresh independent uniform.

Set $Y_i=X_i^2$, $\overline{Y}=n^{-1}{\sum }_i{Y}_i$ and ${\tau }^2={\mathrm{Var}}_F(Y_1)={\mathbb{E}}_F{X}_1^4-{\theta }^2$. The fourth-moment assumption makes ${\tau }^2$ finite, so the variance of a sample mean is $v={\tau }^2/n$. Assume $n\ge 2$, as required for leave-one-out estimation.

The Jackknife variance estimator recomputes the original statistic on each sample with one observation removed. HereBy the definition of the Jackknife variance estimator,For completeness, expanding the centered sum givesConsequently $\mathbb{E}{\hat{v}}_{\mathrm{JACK}}={\tau }^2/n=v$. The jackknife variance estimator is unbiased for this statistic. This calculation uses only the finite fourth moment, not normality of the original observations.

Given the observed data, put $y_i=x_i^2$ and $\overline{y}=n^{-1}{\sum }_i{y}_i$. The empirical distribution assigns mass $1/n$ to each observed $x_i$, including repeated values with their multiplicities. A bootstrap sample consists of $n$ independent draws from this distribution, with replacement. The nonparametric bootstrap variance estimate is the conditional varianceIndependence of the bootstrap draws gives the exact conditional bootstrap variance of a sample mean:Thus its finite-sample expectation is $(n-1)v/n$, in contrast to the unbiased jackknife estimate. This is the jackknife and bootstrap variance of a sample second moment scaling difference.

To approximate it by Monte Carlo simulation, independently generate $B$ bootstrap samples. For replicate $b$, choose $n$ independent uniform indices $I_{b1},\dots ,I_{bn}\in \{1,\dots ,n\}$ and compute $T_b^{\ast }=n^{-1}{\sum }_j{x}_{I_{bj}}^2$. With ${\overline{T}}^{\ast }=B^{-1}{\sum }_b{T}_b^{\ast }$, useFor $B>1$ this is conditionally unbiased for the bootstrap variance and converges to it as $B\to \infty$. Uniform indices can be obtained by $I=\lceil nU\rceil$ for uniforms on $(0,1]$.

Let $G_{\ast }(t)=P_{\ast }(T^{\ast }\le t)$ be the conditional bootstrap distribution of the statistic. A nominal level $1-\alpha$ percentile bootstrap confidence interval isReuse the bootstrap replicates, sort them as $T_{(1)}^{\ast }\le \cdots \le T_{(B)}^{\ast }$, and take the order statistics with ranks $\lceil B\alpha /2\rceil$ and $\lceil B(1-\alpha /2)\rceil$ as empirical quantiles. These are quantiles of the statistic itself, rather than reflected quantiles of its centered error. Increasing $B$ reduces simulation error; it does not turn an approximate confidence interval into an exact finite-sample one.

An analytic alternative follows from the central limit theorem. If ${\tau }^2>0$, then $\sqrt{n}({\hat{\theta }}^n-\theta )\Rightarrow N(0,{\tau }^2)$, and the sample variance $s_Y^2=(n-1{)}^{-1}{\sum }_i(x_i^2-\overline{y}{)}^2$ converges to ${\tau }^2$ by the strong law of large numbers. The Slutsky theorem therefore gives the large-sample normal approximation intervalwhere $z_{1-\alpha /2}$ is the standard normal quantile. Its lower endpoint may be truncated to zero because the parameter is nonnegative. If ${\tau }^2=0$, $X_i^2=\theta$ almost surely and the statistic identifies the parameter exactly; the interval degenerates to that value. The nondegenerate normal interval is asymptotic, not an exact Student interval under arbitrary $F$.

Both the Metropolis–Hastings algorithm and the Gibbs sampler construct a Markov chain with a specified target distribution, usually a posterior distribution known only up to normalization. For a proposal density $q(y\mid x)$ and target density $\pi$, the Metropolis–Hastings acceptance probability isAfter drawing $y$, move there with this probability; otherwise stay at $x$. The unknown normalizing constant of $\pi$ cancels. The original Metropolis algorithm uses a symmetric proposal, $q(y\mid x)=q(x\mid y)$, so the proposal-density ratio disappears. The general Metropolis–Hastings algorithm permits nonsymmetric proposals and needs that correction.

A Gibbs sampler instead chooses a coordinate or block and draws directly from its full conditional distribution, leaving the other coordinates fixed. It is exactly a conditional Metropolis–Hastings update whose proposal is that full conditional. For a change from $x_i$ to $y_i$ with $x_{-i}$ fixed, the acceptance ratio isThus a Gibbs update is an always-accepted Metropolis–Hastings update with the full conditional as proposal. This statement concerns each update: a systematic sweep of several Gibbs kernels preserves the target but need not itself be reversible.

When the full conditional distributions are easy to sample, a Gibbs sampler avoids proposal tuning and rejection, and is particularly convenient for conjugate priors and data augmentation. When a conditional cannot be sampled directly, a Metropolis–Hastings algorithm can use only its evaluable density; such an update can also be inserted into a Gibbs sweep. Block or joint Metropolis–Hastings proposals can move efficiently along strongly correlated directions for which coordinate-wise Gibbs updates mix slowly. Conversely, a badly chosen Metropolis–Hastings proposal can reject most moves or make only tiny moves. Acceptance probability one by itself does not guarantee rapid exploration. Both methods require appropriate irreducibility and convergence checks before their long-run averages are useful posterior estimates.

Write $S_m={\sum }_{i=1}^m{x}_i$ and $S=S_N$. Use shape-rate notation for the Gamma distribution, as in the authoritative PDF: shape $a$ and rate $b$ give a density proportional to $z^{a-1}{e}^{-bz}$. Conditional independence of the Poisson distributions gives a likelihood function proportional toThe product of factorial terms depends only on the data and can be absorbed into the normalizing constant. Multiplying by the two independent Gamma distribution prior kernels and the uniform prior on the possible indices gives the Poisson count change-point model with gamma priors:Because the hyperparameters are positive, each possible index gives finite positive gamma integrals, and summing over the finitely many indices yields a proper posterior distribution. At $m=N$ the second likelihood factor is one, so the posterior law of its unused rate $\varphi$ remains its prior $\mathrm{Gamma}(\gamma ,\delta )$.

The posterior kernel gives two gamma full conditional distributions directly:with rates as the second parameters. For the index, terms independent of the proposed value $j$ cancel. Its conditional weights areDefine ${\ell }_j=j(\varphi -\lambda )+S_j\log (\lambda /\varphi )$ and ${\ell }_{\max }={\max }_j{\ell }_j$. The normalized categorical probabilities areSubtracting the largest log weight prevents needless numerical overflow, an application of the log-sum-exp function.

Precompute the cumulative counts $S_j$. Initialize any $m\in \{1,\dots ,N\}$ and positive rates. At each iteration, draw $\lambda$ from its gamma conditional given the current index, draw $\varphi$ from its gamma conditional given that index, and draw a new index from the categorical probabilities using the newly drawn rates. Repeat and retain the parameter triples after initialization effects have diminished. These are the Gibbs updates for a Poisson count change point.

To verify invariance rather than just identify a sampler, a coordinate update integrates the old joint density as $\pi (\lambda ,\varphi ,m\mid x)=\pi (\lambda \mid \varphi ,m,x)\pi (\varphi ,m\mid x)$ and replaces the first factor by an independent draw from the same conditional. The joint posterior distribution is therefore unchanged. The same argument applies to the other two updates, so their composition preserves the posterior. All gamma densities and all categorical weights are positive on their supports. The full-sweep transition density is consequently positive across the joint support, giving the usual irreducibility and aperiodicity conditions for this proper-posterior Gibbs sampler. This supports posterior long-run averages under the ergodic theorem for a positive Harris recurrent Markov chain.

The event (m=N) represents no change within the observed sequence: every observation belongs to the first-rate segment, and the second-rate segment is empty. It is not the event $\lambda =\varphi$; under continuous independent rate priors that exact equality has posterior probability zero. The endpoint is instead a distinct index value with positive prior mass.

For retained Markov chain Monte Carlo draws $({\lambda }^{(b)},{\varphi }^{(b)},m^{(b)})$, the direct estimate isIts consistency follows from the ergodic theorem for a positive Harris recurrent Markov chain. The draws are correlated, so an independent-binomial standard error would generally be incorrect; uncertainty should account for the Markov chain Monte Carlo asymptotic variance.

A conditional Monte Carlo alternative averages the conditional probability from the preceding solution:The tower property of conditional expectation shows that its posterior expectation is the same target. This replaces an indicator by its conditional expectation; for correlated chains one should assess efficiency rather than automatically assume its entire asymptotic variance is smaller.

An exact finite-sum check is available by integrating out both rates. The marginal posterior weights of a Poisson count change point areThese follow from the Gamma function integral ${\int }_0^{\infty }{z}^{a-1}{e}^{-bz}dz=\Gamma (a)/b^a$. All omitted prior-normalization and factorial factors are independent of $j$, so they cancel from this ratio. This exact calculation provides a useful independent check on either simulation estimate.

Both data augmentation and the expectation-maximization algorithm introduce latent variables $Z$ that make the complete-data model easier to handle. The observed-data likelihood is obtained by summing or integrating $p(x,z\mid \theta )$ over the missing data, and both methods use the conditional law of $Z$ given the observations and current parameters.

In data augmentation for Bayesian inference, one alternates random drawsThese are Gibbs sampler updates on the joint posterior distribution; marginal parameter draws explore the whole posterior, including uncertainty, provided the chain converges appropriately. The parameter conditional incorporates its prior.

The expectation-maximization algorithm is instead a deterministic optimization procedure. Its E step formsand its M step maximizes this function over $\theta$. Exact steps do not decrease the observed-data likelihood function. Their aim is a likelihood stationary point, often a local maximum, rather than draws from a parameter distribution. Adding a log prior to the maximization gives a posterior-mode version, which still does not sample the posterior.

Data augmentation samples conditional missing data and parameters; EM averages the complete-data log likelihood and optimizes parameters. The E step is an expectation of the log likelihood, not the log of an averaged likelihood, and generally is not equivalent to merely replacing each missing value by its conditional mean. A Monte Carlo E-step approximation is a further variant of EM, not automatically the same procedure as a data augmentation chain.

Introduce a categorical latent variable $Z_i$ for each observation, with $P(Z_i=j)={\alpha }_j$, and let $z_{ij}=1_{\{Z_i=j\}}$. Conditional on $Z_i=j$, $X_i$ has a normal distribution of mean ${\mu }_j$ and the shared variance ${\sigma }^2$. The complete-data log-likelihood function isInitialize finite means, positive weights summing to one, and a positive common variance. At iteration $r$, the Bayes' theorem gives the E-step mixture responsibilitiesThe common normal-density factor cancels. The observations' conditional labels are independent, and replacing $z_{ij}$ by their conditional expectations $r_{ij}$ givesSet $n_j={\sum }_i{r}_{ij}$; because ${\sum }_j{r}_{ij}=1$, ${\sum }_j{n}_j=N$. The constrained weight maximization uses a Lagrange multiplier $\eta$: differentiating ${\sum }_j{n}_j\log {\alpha }_j-\eta ({\sum }_j{\alpha }_j-1)$ gives ${\alpha }_j=n_j/\eta$, and summing fixes $\eta =N$. Differentiating the residual quadratic in each mean gives ${\sum }_i{r}_{ij}(x_i-{\mu }_j)=0$. Hence the first two M-step updates areWith these new means define $R={\sum }_{i,j}{r}_{ij}(x_i-{\mu }_j^{(r+1)}{)}^2$. The remaining objective in $u={\sigma }^2$ is $-N\log u/2-R/(2u)$, whose derivative is $-N/(2u)+R/(2u^2)$. For $R>0$ it increases up to $u=R/N$ and decreases afterwards. Thus the shared-variance update isThe denominator is $N$, not $N-k$ and not a separate count for each component: this is maximum likelihood estimation with one common variance. These are the EM for Gaussian mixtures with a common variance updates. Recompute responsibilities with the new parameters and repeat until the observed-data likelihood changes sufficiently little.

The monotonicity claim also follows directly. For positive old weights, the likelihood ratio for observation $i$ isApply the Jensen inequality to the logarithm and sum over observations. This yieldsfor the maximizing M step. This proves likelihood increase, not convergence to the global maximum. Use several initializations to reduce sensitivity to local maxima. Compute responsibilities with a log-sum-exp function normalization to avoid numerical underflow.

If $n_j=0$, that component's updated weight is zero and its mean is unidentified; retain its old mean or remove the empty component. Strictly positive interior initialization avoids this case in exact arithmetic while the variance stays positive. If $R=0$, there is no positive finite variance maximizing this M step: the objective tends to infinity as ${\sigma }^2\downarrow 0$. Such boundary degeneracy needs a variance constraint or regularization, as discussed in the next solution. For $k=1$ the formulas reduce to the ordinary sample mean and the maximum-likelihood variance $N^{-1}{\sum }_i(x_i-\overline{x}{)}^2$.

Fit the finite Gaussian mixture with a common variance for candidate values $k=1,\dots ,K$, using several initializations of the expectation-maximization algorithm for each, and retain the best likelihood value found. Maximized likelihood alone cannot select the number of components: a smaller mixture embeds in a larger one by giving an added component zero weight, so its likelihood cannot decrease.

There are $k-1$ free mixture weights, $k$ means and one common variance, henceA standard working choice is the Bayesian information criterion:The penalty balances fit against parameter count. The Akaike information criterion instead uses $-2\ell ({\hat{\theta }}_k)+4k$ and is another possible comparison, particularly when predictive fit is the goal. The parameter count would differ if component variances were separate, but that is not this model.

Likelihood fitting must first be made well posed. Even a shared variance does not prevent all degeneracy: if $k=N$, choose ${\mu }_j=x_j$, ${\alpha }_j=1/N$ and let the common variance tend to zero. Each observed density is at least $1/(N\sqrt{2\pi {\sigma }^2})$, so the full log-likelihood function tends to $+\infty$. The same construction works when $k$ reaches the number of distinct observed values. Thus one must not feed infinite fitted likelihoods into an information criterion. Restrict candidate sizes to a sensible range below that degeneracy, or impose a fixed positive variance floor and use that same constrained model in every comparison. Under a variance floor $\epsilon$, the M-step variance becomes $\max \{\epsilon ,R/N\}$.

There is also nonregular mixture model selection: at a missing component its weight is on a boundary and its mean is unidentified, while permuting labels leaves the distribution unchanged. Consequently an ordinary chi-squared calibration for a likelihood-ratio test between successive component counts is not generally justified. Information criteria are practical comparisons, not an automatic universal proof that the true number has been recovered. Predictive cross-validation or a suitably constructed parametric bootstrap comparison can supplement them. Select the component count by a penalized, well-defined mixture fit, rather than by unpenalized likelihood increase.