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

Put ${\theta }_0=1$ and define the autoregressive polynomial and moving-average polynomial byFor a minimal, coprime ARMA representation, the usual simultaneous causality and invertibility root criteria for an ARMA model areThe autoregressive condition makes $\Theta /\Phi$ a stable one-sided filter and gives the causal second-order stationary process. The moving-average condition makes its inverse a stable one-sided filter, so the white noise can be reconstructed from present and past observations. If the two polynomials have a common factor, first use common factor cancellation in an ARMA model; necessary-and-sufficient root statements apply to the reduced noise-driven representation. Stationarity here is second-order stationarity; weak white noise alone does not assert strict stationarity of all finite-dimensional laws.

For $p=1$, write $\varphi ={\varphi }_1$, so $|\varphi |<1$. Expanding $(1-\varphi z{)}^{-1}={\sum }_{r\ge 0}{\varphi }^r{z}^r$ gives the Wold representationIn particular $c_0=1$, $c_j=\varphi c_{j-1}+{\theta }_j$ for $1\le j\le q$, and $c_j={\varphi }^{j-q}{c}_q$ for $j\ge q$. The coefficients are square-summable, so this series converges in mean square. Causality and invertibility identify the driver with the linear innovation process, justifying the Wold terminology. For $\varphi =0$, powers with exponent zero equal one and the formula reduces to the finite moving average.

To obtain the covariance equations, subtract $\varphi X_{t-1}$ and take covariance with $X_{t-k}$. Orthogonality of the white noise and the causal expansion implyThus the innovation coefficients of an ARMA(1,q) process give the entire triangular system in one formula:At $k=0$, use ${\gamma }_{-1}={\gamma }_1$, since the autocovariance of a real stationary process is symmetric. The terms are consequently ${\sigma }^2(c_0+{\theta }_1{c}_1+\cdots +{\theta }_q{c}_q)$ at zero, ${\sigma }^2{\theta }_q{c}_0$ at $q$, and the intervening partial sums as required.

For $q=1$, put $\theta ={\theta }_1$ and $c_1=\varphi +\theta$. The equations becomeSubstitution of the second into the first, followed by the geometric recursion, yields the autocovariance of a causal ARMA(1,1) process:This includes $\varphi =0$: only lags $\pm 1$ survive, as for an MA(1) process.

The state-space model as printed supplies only uncorrelated white noise. That is insufficient for the claimed normal conditional distributions. For example, take $S_0$ normal, independent process noises $w_t$ taking values $\pm \sqrt{W}$ equally often, and independent Gaussian observation noises. With $G=1$ and no previous observations, $S_1=S_0+w_1$ is a two-component shifted mixture distribution, not normal for $W>0$. Its characteristic function is $e^{i\mu s-P_0{s}^2/2}\cos (\sqrt{W}s)$, which has zeros, whereas a normal characteristic function has none. Similarly a non-Gaussian observation noise prevents the asserted conditional normal law for $X_t$.

The intended scalar Gaussian Kalman recursion requires $w_t\sim N(0,W)$ and $v_t\sim N(0,V)$ independent of each other and of the previous state and observation history. Assume these Gaussian independence conditions, $V>0$, $W\ge 0$, and the given normal filtering law. Conditional on the past, $G{S}_{t-1}+w_t$ is a sum of independent normals, henceThe new observation noise is independent of $S_t$ and the past, so $X_t\mid S_t,{\mathcal{F}}_{t-1}\sim N(F{S}_t,V)$. Write $S_t=m_t+{\eta }_t$ and $X_t=F{m}_t+F{\eta }_t+v_t$. The independent normal pair $({\eta }_t,v_t)$ shows directly thatThe covariance $F{R}_t$ comes from the common ${\eta }_t$ term. Applying the conditional multivariate normal distribution formula with $X_t$ as the conditioned coordinate gives the Kalman filter updateThe final expression verifies nonnegativity of the posterior variance.

For $F=G=1$, defineThe recursions become ${\hat{S}}_t=(1-{\alpha }_t){\hat{S}}_{t-1}+{\alpha }_t{X}_t$ and $P_t={\alpha }_{t}V$. If $P_t\to P$, then ${\alpha }_t\to \alpha =P/V$ and, with $c=W/V$,The nonnegative solution is the steady-state local-level Kalman gain:For $c=0$ the limit is zero; for $c>0$ it lies strictly between zero and one. These are exact conditional-law results under the additional Gaussian assumptions, not consequences of the weaker printed noise description alone.

Let $h\ge 0$ be proportional to the target probability density function, with $0<Z=\int h(x)dx<\infty$. Define the ratio-of-uniforms method through the regionDraw $(U,V)$ uniformly with respect to area on $A_h$ and return $X=V/U$. To prove its law, change variables from $(u,x)$ to $(u,v)=(u,ux)$. The Jacobian determinant is $u$, and the allowed interval for $u$ is $0<u\le \sqrt{h(x)}$. ConsequentlyThe joint density of $(U,X)$ is therefore $2u/Z$ on that interval. Integrating out $u$ givesFor an already normalized density, $Z=1$ and the region has area $1/2$. The name refers to the two coordinates of a uniformly sampled region; drawing two unrelated uniforms on arbitrary intervals and taking their ratio would not give the target law.

The region $A_h$ usually has a curved boundary, so sampling uniformly inside it directly is inconvenient. Rejection sampling replaces it by an easily sampled envelope, usually a rectangle. Choose finite bounds $a,b_{-},b_{+}$ withThen $A_h\subseteq (0,a)\times (b_{-},b_{+})$ up to null boundaries. A general rectangle algorithm uses only independent uniform deviates: draw $U_1,U_2$ on $(0,1)$, set $u=a{U}_1$, $v=b_{-}+(b_{+}-b_{-})U_2$, and return $v/u$ if $u^2\le h(v/u)$; otherwise repeat. If a numerical generator returns $u=0$, discard that draw before division. No knowledge of $Z$ is required to run the method. Its acceptance rate is the ratio of the two areas:More generally the same area ratio applies to any finite-area envelope $B$ sampled uniformly. A change in the proportionality constant of $h$ scales both ratio coordinates and their areas consistently.

Finite rectangle bounds are a real condition, not a consequence of integrability alone. For instance $h(x)=(1+|x|{)}^{-3/2}$ is integrable, but $|x|\sqrt{h(x)}$ is unbounded. A finite-area envelope for a ratio-of-uniforms region gives a general extension: choose a finite-area open envelope containing $A_h$, decompose it up to null boundaries into countably many disjoint rectangles, and select a rectangle with probability proportional to its area using a uniform draw and cumulative probabilities. Use two more uniforms to sample within it, then apply the same acceptance test. This samples the envelope uniformly and has acceptance $Z/(2\mathrm{area}(B))$. Such envelopes exist for every finite-area measurable $A_h$; obtaining usable bounds or an effective rectangle decomposition is density-specific. The finite-rectangle version above is the usual practical implementation.

For $\sigma >0$, use the unnormalized normal density $h(x)=e^{-x^2/(2{\sigma }^2)}$. Its ratio region isThe bound on $u$ is one because $h$ is maximized at zero. To bound $|v|$, maximize $x\sqrt{h(x)}=x{e}^{-x^2/(4{\sigma }^2)}$ for $x\ge 0$. Its derivative vanishes at $x=\sqrt{2}\sigma$, giving $b=\sigma \sqrt{2/e}$. Thus the normal ratio-of-uniforms envelope is $0<u<1$, $-b<v<b$.

Draw independent $U_1,U_2\sim \mathrm{Unif}(0,1)$, set $u=U_1$ and $v=b(2U_2-1)$, and accept precisely whenreturning $X=v/u$. Equivalently compare logarithms to avoid underflow. Since $Z=\sigma \sqrt{2\pi }$, the region area is $\sigma \sqrt{\pi /2}$ and the rectangle area is $2\sigma \sqrt{2/e}$. ThereforeIt is independent of the normal scale. In the diagram the vertical coordinate is normalized by $\sigma$, so the same shape works for every positive scale.

For independent standard Cauchy distribution draws $X_1,\dots ,X_n$, the code averages the indicators $I_i=1_{\{X_i>2\}}$. Each is a Bernoulli random variable with parameterHence it is an unbiased estimator of that tail probability, and independence givesThe variance is finite even though the underlying Cauchy draws have no finite mean: the variables being averaged are bounded indicators, not the Cauchy values themselves.

For $Y=2/(1-U)$ with $U$ uniform on $(0,1)$,Thus $Y$ has a Pareto distribution and is the importance sampling proposal, concentrated on the required tail. The target tail integrand is $f(y)=1/(\pi (1+y^2))$, so its importance weight isThe function returns $\tilde{\theta }=n^{-1}{\sum }_{i}w(Y_i)$. Since ${\mathbb{E}}_{q}w(Y)={\int }_2^{\infty }f(y)dy=\theta$, it is an unbiased estimator. For its second moment,An antiderivative is $\frac{1}{2}(\arctan y-y/(1+y^2))$, so the integral is $\frac{1}{2}(\pi \theta +2/5)$. Independence of the proposals givesThis is importance sampling of a Cauchy tail, with weights between $2/(5\pi )$ and $1/(2\pi )$. Their narrow range explains the much smaller variance than the indicator estimate, without appealing to nonexistent Cauchy moments.

There is a genuine normalization error in the printed function. Write $Y_i\sim \mathrm{Unif}(0,2)$ independently, $g(y)=1/(\pi (1+y^2))$ and $H(y)=y^2-4/3$. The internal sum produces a scalar; taking its mean does not divide by $n$. The literal returned value is thereforeSince $\mathbb{E}H=0$ and $\mathbb{E}g=\arctan (2)/(2\pi )=1/4-\theta /2$, its expectation isIt does not estimate $\theta$ unbiasedly unless $n=1$, and no fixed control variate coefficient removes this bias. In particular the actual call with $\beta =0$ returnsThis is an instance of summing before taking a mean changes Monte Carlo normalization.

The intended corrected method replaces the sum by a sample mean:It uses symmetry to write the tail as $\theta =1/2-{\int }_0^{2}g(y)dy$, then Monte Carlo integration on a bounded interval, and finally a mean-zero control variate. For $\beta =0$ the intended estimator is $1/2-(2/n){\sum }_{i}g(Y_i)$, whose summands are smooth and bounded. Keeping the average normalization, $\mathbb{E}{T}_{\beta }=\theta$ andDifferentiating this quadratic gives the Cauchy tail control variate on a bounded interval coefficientHere $\mathbb{E}{Y}^2=4/3$, $\mathbb{E}{Y}^4=16/5$, so $\mathrm{Var}(H)=64/45$. AlsoConsequentlyThe negative sign is consistent with $g$ decreasing and $H$ increasing. One may estimate the covariance and variance using a separate pilot sample and fix the resulting coefficient for a fresh sample. This coefficient minimizes the literal function's variance as well, but its missing factor of $1/n$ still leaves the bias above.

For a target probability density function $\pi$ and a proposal distribution with density $q(x^{\prime }\mid x)$, the Metropolis–Hastings algorithm starts in the target support and repeats: propose $x^{\prime }\sim q(\cdot \mid x)$, draw an independent uniform $U$, and move to $x^{\prime }$ ifOtherwise retain $x$. A zero reverse-proposal density makes the acceptance zero. Unknown multiplicative constants of $\pi$ cancel. For distinct states,which is symmetric in the states. Together with the probability of holding at the current state, this proves detailed balance and invariance of $\pi$. Successive values form a dependent Markov chain; convergence from arbitrary initial states additionally needs the appropriate irreducibility and aperiodicity, not just this invariant-density identity.

The Gibbs sampler updates a coordinate using its full conditional distribution while keeping the others fixed. For coordinate $i$, write $r=x_{-i}$ and $\pi (x)={\pi }_{-i}(r){\pi }_i(x_i\mid r)$. Its proposal replaces $x_i$ by $x_i^{\prime }\sim {\pi }_i(\cdot \mid r)$. The MH ratio, interpreted on this one-coordinate proposal space, isThus every Gibbs coordinate update is an MH update with acceptance one. A random choice of coordinate gives a mixture of these invariant kernels; a deterministic sweep gives their composition. Both preserve $\pi$. The latter composition need not be reversible, so the special-case statement is about the coordinate kernels, not a claim that every entire systematic sweep is one reversible MH kernel.

At $\sigma =1$, complete the square:For fixed $y$, all factors not involving $x$ enter its conditional normalizing constant. Interchanging the roles of the coordinates givesThese are the full conditional distributions of the standard bivariate normal distribution with correlation $\rho$.

A systematic Gibbs sampler starts at a chosen $(x_0,y_0)$ and, at each sweep, draws independent standard normals $Z_{1,n},Z_{2,n}$ and setsThe second update uses the newly drawn $x_{n+1}$. Both coordinate updates preserve the target by part (a), hence so does the sweep. The Gaussian two-coordinate Gibbs recursion is explicitlyan AR(1) recursion with innovation variance $1-{\rho }^4$. Its stationary variance is one, and its correlation at sweep lag $k$ is ${\rho }^{2k}$, displaying the dependence of the sample.

The printed general-$\sigma$ expression has a separate normalization defect: for $\sigma >0$ its integral is $\sigma$, not one, since the exponent has covariance ${\sigma }^2\left(\begin{array}{cc}1& \rho \\ \rho & 1\end{array}\right)$, whose normalizing denominator is $2\pi {\sigma }^2\sqrt{1-{\rho }^2}$. Retaining that exponent requires ${\sigma }^2$ in the prefactor. The requested $\sigma =1$ target is already normalized, so this defect does not alter the conditional distributions above.

Use the random-walk proposal specified by its centre and covariance: for the current state $z=(x,y)$ draw $Z_1,Z_2$ independently standard normal and setIts properly normalized density of the proposal distribution isIt is symmetric: $q(z^{\prime }\mid z)=q(z\mid z^{\prime })$. This is a centred-at-the-current-state proposal, rather than an unshifted draw from the density notation in the source. The random-walk Metropolis algorithm consequently needs only the target ratio. Put $H(x,y)=x^2+y^2-2\rho xy$. ThenAn independent uniform decides the move; on rejection both coordinates remain at their previous values. A proposal with lower $H$ is always accepted. This construction uses the joint density directly, without needing either conditional distribution. Its continuous positive proposal density supplies communication over the target support, and the holding probability avoids a deterministic periodic transition pattern.

Tune the proposal scale using pilot runs on a grid of candidate ${\sigma }_q$ values. After allowing for initial transients, compare the effective sample size of a Markov chain per unit computation time for representative functions, such as both coordinates and a quadratic function. Equivalently inspect their integrated autocorrelation time, using $\mathrm{ESS}\approx n/(1+2{\sum }_{k\ge 1}\mathrm{Corr}(h(Z_0),h(Z_k)))$. Choose a useful scale and keep it fixed for the production run. Another pilot criterion is the average squared accepted displacement, possibly in coordinates standardized by the target covariance, which balances move size and rejection.

Acceptance alone is not an efficiency criterion. As ${\sigma }_q\downarrow 0$, proposals are almost the current state, so acceptance tends to one while exploration becomes arbitrarily slow. Very large scales propose mostly low-density states and produce long sequences of identical values. The proposal scale and random-walk Metropolis efficiency trade-off requires both movement and acceptance; neither maximizing nor merely reporting the acceptance probability guarantees a useful sample. With strong correlation, an isotropic proposal also struggles with the target's narrow direction; covariance-scaled proposals can improve this geometry beyond tuning one scalar variance.

The observed-data likelihood function integrates out the latent variable:Choose an initial parameter ${\theta }^{(0)}$ with positive likelihood. At iteration $t$, the E-step of the EM algorithm forms the entire expected complete-data log-likelihoodThe expectation is under the conditional missing-data law at the old parameter, while the argument $\theta$ remains free. The M-step then choosesThe observed data $x$ stay fixed throughout. This is not obtained merely by substituting the conditional mean of $z$ into a nonlinear log-likelihood: the expectation of all required complete-data terms must be taken. A generalized EM step only requires $Q({\theta }^{(t+1)}\mid {\theta }^{(t)})\ge Q({\theta }^{(t)}\mid {\theta }^{(t)})$, which suffices for the monotonicity proved next. Exact maximization assumes that an admissible maximizer exists.

Write ${\theta }_0={\theta }^{(t)}$ and $h_0(z)=f(x,z;{\theta }_0)/L(x\mid {\theta }_0)$. Assume the likelihoods are positive and the expected log likelihoods used in the step are finite. On the support where $h_0>0$, let $r(z)=f(x,z;\theta )/f(x,z;{\theta }_0)$. Integrating over that support givesIt is an equality when the candidate joint density has no additional mass outside the old support. Concavity of the logarithm and the Jensen inequality now implyThus the argument does not silently require an identical support at every parameter. The M-step selects a value whose right-hand side is nonnegative, proving EM likelihood monotonicity:When the conditional densities allow finite divergence, the same calculation has the informative exact formThe Kullback-Leibler divergence is nonnegative, since ${\mathbb{E}}_{h_0}\log (h_{\theta }/h_0)\le \log {\int }_{\{h_0>0\}}{h}_{\theta }\le 0$. This also explains why an increase of the surrogate objective guarantees an increase of the observed objective, while equality of the two likelihoods is possible.

EM likelihood monotonicity makes EM an ascent method, not a global-optimization theorem. If the log-likelihood is bounded above along the iterates, its nondecreasing values converge. That alone does not show that the parameters converge, nor that any limiting parameter is a global maximum-likelihood estimator. Under additional regularity, accumulation points can be stationary; local maxima, other stationary points or flat sets can depend on initialization. Multiple starting values and comparison of the achieved observed likelihoods are therefore useful.

The upper-bound hypothesis is substantive. In a finite mixture model with a freely shrinking normal-component variance, placing one component at an observed value can make the likelihood unbounded while other components cover the remaining observations. In that situation monotonicity provides no finite maximum to converge to. EM may also improve very slowly, and a small change in likelihood alone is not proof of an accurate global solution. The guarantee in part (b) is nondecrease, not strict increase at every iteration.

Both EM algorithm and data augmentation exploit a latent completion $z$, but they use it differently. EM takes a conditional expectation of the complete-data log-likelihood and then maximizes over the parameter, producing a deterministic sequence of point estimates when the M-step is uniquely specified.

In Bayesian data augmentation, choose a prior distribution $p(\theta )$ and targetAlternate the Gibbs sampler updatesThese are stochastic draws from full conditional distributions, not an expectation followed by maximization. Their invariant law is the augmented joint posterior, and discarding $z$ gives posterior inference for $\theta$, including uncertainty rather than only an MLE. A sampled completion followed by maximization would instead be a stochastic EM variant. Individual likelihood values along a data-augmentation chain need not increase: the chain is designed to sample a posterior, not to ascend the likelihood.