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

For a weakly stationary process with mean $m$ and nonzero variance, put $\gamma (h)=\mathrm{Cov}(X_{t+h},X_t)$ and $\rho (h)=\gamma (h)/\gamma (0)$. The function $\rho$ is its autocorrelation function. For observations $x_1,\dots ,x_N$, one common sample autocorrelation function isA correlogram plots these estimates against lag; the population analogue plots $\rho (h)$.

The lag-$h$ partial autocorrelation function measures the correlation between $X_t$ and $X_{t-h}$ after removing their linear projections on the intervening observations. Its sample counterpart can be calculated by the Yule-Walker equations: solveand take $\hat{\alpha }(h)={\hat{\varphi }}_{h,h}$, assuming the matrix is nonsingular. This defines a standard sample partial autocorrelation function.

For a minimal causal autoregressive model of order $p$, the autocorrelation typically decays geometrically, possibly with damped oscillations, whereas the partial autocorrelation function is zero beyond lag $p$: after those $p$ preceding values have been projected out, the innovation is orthogonal to older values. For a moving-average model of order $q$, nonoverlapping sets of driving white noise give $\gamma (h)=0$ for $|h|>q$; its partial autocorrelation function generally tails off when the representation is invertible. Thus an ACF cutoff suggests an MA order; a PACF cutoff suggests an AR order. An ARMA model generally has neither cutoff. These are population properties; sample noise makes cutoffs approximate, so fitted residuals and uncertainty should also be inspected. Under independent white noise, fixed-lag sample correlations have approximate standard error $N^{-1/2}$, giving a rough diagnostic band rather than a universal simultaneous test.

For a discrete-time weakly stationary process, the spectral density of a stationary process, when it exists, is a nonnegative function satisfyingIf ${\sum }_h|\gamma (h)|<\infty$, thenA general second-order stationary process has a spectral measure of a stationary time series, which need not have a spectral density: a centered random constant has an atom at zero. Thus the spectral density calculation below assumes that $f_X$ exists; the corresponding measure identity holds without that extra assumption.

Absolute summability of the filter coefficients gives mean-square convergence becauseConsequently the bilateral filtered process is well defined, with constant mean $m{\sum }_s{a}_s$. Its covariance isThe sum is absolutely convergent, since $|{\gamma }_X(j)|\le {\gamma }_X(0)$ and $\sum |a_s|<\infty$. It depends only on lag, proving weak stationarity of $Y$.

Insert the spectral representation and interchange summation and integration using absolute summability and the finite total spectral mass. For real coefficients,Therefore the spectral density transformation under a linear filter isFor a bilateral sequence, $A$ is guaranteed to converge on the unit circle. Its negative powers need not converge inside the disk: for example, $a_{-j}=2^{-j}$ makes the series diverge at $z=1/4$. Only the unit-circle values are used here, so the unnecessarily broad domain printed for $A$ does not affect the proof.

Center the five coefficients at lags $-2,-1,0,1,2$. Their transfer function is real:Choosing consecutive causal lags instead only adds a phase and leaves the squared modulus unchanged. After $k$ applications, the white noise spectral density has becomeWrite $c=\cos \omega$. Completing the square givesFor $-1\le c\le 1$, its range is $[-1/3,7/6]$, and its absolute maximum is $7/6$, attained exactly when $c=1/2$. Consequently, on the conventional nonnegative frequency range $[0,\pi ]$,On the two-sided range $[-\pi ,\pi ]$ there is also a peak at $-\pi /3$; the ratio there is one for every $k$. This is a counterexample to the claim if read literally on a signed frequency domain. The correct two-sided exclusion is $\omega \ne \pm \pi /3$.

The result illustrates how repeated symmetric filtering can amplify an oscillation. Although $H(0)=1$ preserves a constant component, $H(\pi /3)=7/6>1$ amplifies period-six oscillations. Repetition selects an increasingly narrow frequency band around those peaks; it does not simply smooth all noise away. In fact the unnormalized output variance eventually grows without bound, since a fixed interval around either peak has $|H|>1$.

Difference the local-level state-space model to eliminate its latent level:Independence and Gaussianity give a stationary Gaussian process withFor $V,W>0$, letThese satisfy $q(1+{\vartheta }^2)=W+2V$ and $q\vartheta =-V$. Therefore $D$ has the covariance of the moving-average model ${\epsilon }_t+\vartheta {\epsilon }_{t-1}$, where $\epsilon$ is Gaussian white noise of variance $q$. This is an actual representation, not just covariance matching: define ${\epsilon }_t={\sum }_{j\ge 0}(-\vartheta {)}^j{D}_{t-j}$. The convergent inverse filter has constant spectral density $q/(2\pi )$, so its Gaussian coordinates are independent, and multiplication by $1+\vartheta B$ recovers $D$.

ThusThis is the formal ARMA equation with autoregressive coefficient one. Its unit root means that the levels are generally not stationary: if the initial state has finite variance and is independent of future noises, $\mathrm{Var}(S_t)=\mathrm{Var}(S_0)+tW$. In stationary-model terminology the correct classification is ARIMA(0,1,1). The requested ARMA description must therefore allow a unit-root equation.

When $V=0$, the differences are just $w_t$; when $W=0$ and $V>0$, they are $v_t-v_{t-1}$, a noninvertible MA(1) with coefficient $-1$. If both variances vanish, the model is deterministic apart from a possible initial random level. These limits explain the qualifications on the invertible representation above.

Conditionally on past observations, adding the independent state increment gives the predictionThe new observation then has conditional distribution $N({\hat{S}}_{t-1},R_t+V)$, and its conditional covariance with $S_t$ is $R_t$. Equivalently, multiplying the state prediction probability density function by the observation likelihood and completing the square givesFor positive variances, rearrangement yields the scalar Gaussian Kalman recursionHere the innovation is the observation minus its predicted mean. The displayed covariance formula also handles zero variances whenever its denominator is nonzero; a completely deterministic prediction and observation require no probabilistic update. Independence of the new noises from the past and the initial state is the usual state-space assumption needed for this conditional calculation.

A constant posterior variance must be a fixed point of the recursion in part (b). Multiplying $P=V(P+W)/(P+W+V)$ by its positive denominator givesConversely, if a nonnegative $P$ solves this equation and the filter is initialized with $P_0=P$, the recursion preserves $P_t=P$ at every step. The fixed-point equation alone does not force an arbitrary initial covariance to be constant.

For $W,V>0$, the admissible root and its steady-state local-level Kalman gain areThe state estimate then obeys exponential smoothing:Iteration makes the weights explicit:Thus older observations receive geometrically decreasing weights. Consistently with part (a), the innovations variance is $q=P+W+V$ and the differenced MA coefficient is $\vartheta =-(1-K)$. The degeneracies $W=0,P=0$ and $V=0,P=0$ give gain zero and gain one respectively when the denominator is positive.

Treat the supplied uniform variates as independent $U(0,1)$ draws, as required for Monte Carlo simulation. Use $n$ of them and returnEach indicator has a Bernoulli distribution of parameter $p$, and the indicators are independent. For any subset of $r$ successes the probability is $p^r(1-p{)}^{n-r}$; there are $\left(\genfrac{}{}{0ex}{}{n}{r}\right)$ such subsets. Thus $\mathbb{P}(X=r)=\left(\genfrac{}{}{0ex}{}{n}{r}\right)p^r(1-p{)}^{n-r}$, the binomial distribution. The construction includes the cases $p=0,1$ and $n=0$.

For rate $\lambda >0$, returnFor $x\ge 0$, $\mathbb{P}(X>x)=\mathbb{P}(U_1<e^{-\lambda x})=e^{-\lambda x}$, so $X$ has the exponential distribution of rate $\lambda$. This is inverse transform sampling; using $-\log (1-U_1)/\lambda$ is equivalent. Uniform draws at exact endpoints can be excluded, since those events have probability zero under the ideal continuous law.

Use uniform-envelope rejection sampling for a beta distribution. Put $r(x)=x^{a-1}(1-x{)}^{b-1}$ for $0<x<1$ and choose $M=\sup r$. For $a,b>1$, differentiation of $\log r$ gives the mode $(a-1)/(a+b-2)$, soIf $a=1$ or $b=1$, take $M=1$; this also covers $a=b=1$. Draw $X=U_{2j-1}$ and accept it ifOtherwise continue with the next independent pair. The probability of acceptance is ${\int }_0^{1}r(x)dx/M=B(a,b)/M>0$, so the procedure terminates almost surely. Conditional on acceptance, the probability density function is proportional to $r(x)$, hence is exactly the beta distribution of parameters $a,b$. The restriction $a,b\ge 1$ guarantees the bounded uniform envelope.

For $\sigma >0$, use the unnormalized normal shape $h(x)=e^{-x^2/(2{\sigma }^2)}$. The ratio-of-uniforms method samples uniformly fromIts boundary is $|v|\le 2\sigma u\sqrt{-\log u}$. The largest possible $|v|$ occurs at $u=e^{-1/2}$ and equals $b=\sigma \sqrt{2/e}$, giving the normal ratio-of-uniforms envelope $0<u<1$, $-b<v<b$.

Use independent uniform pairs to propose $u=U_{2j-1}$ and $v=b(2U_{2j}-1)$. Accept precisely whenTo verify its distribution, change coordinates from $(u,v)$ to $(u,x)$, with $v=ux$ and Jacobian determinant $u$. The accepted point is uniform on $\mathcal{A}$, and the marginal probability density function of $x$ is proportional toThus the result is $N(0,{\sigma }^2)$. The region has positive finite area $\sigma \sqrt{\pi /2}$, and the acceptance probability is $\sqrt{\pi e}/4$. The case $\sigma =0$ requires simply returning zero.

Assume the draws are independent from $g$, that the target-weighted integrand is covered by its support, and that $\int |\theta (x)|f(x)dx<\infty$. Define the importance weight $w(x)=f(x)/g(x)$. The change-of-density identity isTherefore importance sampling uses the unbiased estimatorThe strong law of large numbers gives consistency. A finite second moment of the weighted integrand additionally gives finite variance and a usual independent-sample central limit theorem. A useful proposal places enough probability where $f|\theta |$ is large and avoids tiny proposal probability density function there; common support alone does not guarantee finite variance.

Write $f=f_0/C_f$ and $g=g_0/C_g$, where $f_0,g_0$ are known nonnegative kernels. If both constants are known, the ordinary weight is $(C_g/C_f)f_0/g_0$. If their ratio is unknown, use self-normalized importance sampling:Indeed, under the normalized proposal law,The strong law of large numbers makes their empirical ratio converge to $\mu$, provided the numerator is absolutely integrable and $0<C_f,C_g<\infty$. Multiplicative constants cancel, whether the unknown one belongs to the target, proposal, or both. Unlike the estimator in part (i), this ratio is generally biased at finite sample size. Exact sampling from $g$ is still presumed; knowing only an unnormalized proposal does not itself supply a sampler.

Set $H_i=f(x_i)\theta (x_i)/g(x_i)$. Independence gives $\mathrm{Var}(n^{-1}{\sum }_i{H}_i)=\mathrm{Var}(H_1)/n$, while part (i) gives ${\mathbb{E}}_g{H}_1=\mu$. Its second moment isConsequentlyThe variance is finite exactly when the displayed second-moment integral is finite, assuming $\mu$ exists. An infinite integral means infinite variance, not failure of the change-of-density identity.

For the standard Cauchy distribution tail, use the indicatorThen $\int \theta (x)f(x)dx$ is the desired probability, and the ordinary importance sampling estimator is $n^{-1}{\sum }_{i}f(x_i)1_{\{x_i\ge k\}}/g(x_i)$. The proposal must cover the nonzero integrand, namely the upper tail. A proposal restricted to $[0,k]$ cannot directly estimate this indicator integral; part (v) instead uses a complementary bounded-interval identity.

Take $k>0$, so the proposal is the uniform distribution on $[0,k]$, with probability density function $g(x)=1/k$. Symmetry of the standard Cauchy distribution givesThus the proposed expression is exactly the unbiased importance sampling estimator of this complementary integral:It has finite variance because $f$ is bounded on $[0,k]$. Explicitly,These follow from the antiderivatives of $(1+x^2{)}^{-1}$ and $(1+x^2{)}^{-2}$. The target itself is $\mu =1/2-\arctan (k)/\pi$. At $k=0$ it is exactly $1/2$ without simulation; for negative $k$ the printed uniform-interval construction must be changed.

Antithetic variates preserve each draw's marginal law but couple two estimates negatively. For an average of identically distributed $H,H^{\prime }$,Negative covariance improves on two independent evaluations at the same computational budget.

Here generate independent $U_1,\dots ,U_r\sim U(0,1)$ and pair $x_j=k{U}_j$ with $x_j^{\prime }=k(1-U_j)$. Both are uniform on $[0,k]$. The antithetic Cauchy-tail integration on a finite interval estimator isIt is unbiased. Put $h(u)=f(ku)$. This decreases, whereas $h(1-u)$ increases. For an independent copy $U^{\prime }$ of $U$, the covariance identity givesEach product is nonpositive, and for $k>0$ it is strictly negative off the diagonal. Thereforewhich is the variance using $2r$ independent evaluations in part (v). Thus the comparison accounts for the doubled number of probability density function evaluations in each pair.

The linear predictor is unchanged by ${\alpha }_j\mapsto {\alpha }_j+c$, ${\beta }_t\mapsto {\beta }_t+d$, $\mu \mapsto \mu -c-d$. Without constraints, different parameter triples therefore give the same rates and likelihood. Identifiability requires choosing one representative of this redundancy.

Taking the first levels as reference yields ${\alpha }_1={\beta }_1=0$. ThenThese identify all three free parameters, while ${\lambda }_{22}={\lambda }_{21}{\lambda }_{12}/{\lambda }_{11}$ is the model's no-interaction constraint. The zero reference effects fix a parametrization; they do not remove the baseline rate or impose equal rates across groups.

Use the standard interpretation that the counts are conditionally independent given the parameters, with mutually independent stated priors. Let $a_2={\alpha }_2$, $b_2={\beta }_2$, and defineUnder the reference constraints, the Poisson distribution likelihood, up to data-only factors, isMultiplying by the shape-rate gamma prior gives the Poisson-gamma conjugacy with unequal exposuresHere $b$ is a rate, matching the probability density function convention on the PDF's first page. The factor $I$ represents the number of replicates per cell, not the total count.

An exact Gibbs sampler update is available despite the non-Gaussian intercept: draw ${\theta }^{\prime }$ from the conditional gamma distribution law in part (ii), then setTogether with conditional updates of the two remaining effects, this targets the required joint posterior in the log-intercept coordinates. The log-gamma prior for a Poisson log-intercept includes the change-of-variable factor: with $d\theta /d\mu =e^{\mu }$,Using $(a+T-1)\mu$ instead would omit the Jacobian and give the wrong conditional law for $\mu$.

The full conditional log probability density functions of the effects are, up to constants,These are Gaussian-prior Poisson log-effect conditional probability density functions. The exponential likelihood term prevents a standard normal distribution conjugate update. A Metropolis-within-Gibbs algorithm conveniently updates one effect at a time without computing its conditional normalizing constant.

For example, propose $a_2^{\prime }=a_2+s_{\alpha }Z$ with $Z\sim N(0,1)$, and accept withUse an analogous normal distribution random-walk proposal for $b_2$. The proposal is symmetric and lives on the appropriate unrestricted real parameter space, so no proposal ratio is needed. It makes local moves compatible with a unimodal conditional probability density function. Tune the positive scale in a pilot or freeze it after warm-up; a scale far too large causes rejection and a scale far too small moves slowly.

A useful scale can be inferred from the conditional curvature: ${\ell }_{\alpha }^{\prime \prime }(a_2)=-I\theta (1+e^{b_2})e^{a_2}-{\sigma }_1^{-2}<0$, with the analogous formula for $b_2$. A normal proposal near the conditional mode with variance approximately $-1/{\ell }^{\prime \prime }$ is another sensible choice, but an asymmetric or state-dependent proposal must include its full reverse-to-forward proposal-density ratio. Metropolis–Hastings is convenient, not logically mandatory; these concave conditionals also admit other specialized samplers.

Introduce a model indicator $M\in \{0,1\}$. In $M=0$, set $\mu =0$, hence $\theta =1$, and retain the two group effects. In $M=1$, include $\mu$ with the prior induced by the gamma baseline. Choose positive prior model probabilities $p_0,p_1$ summing to one; the supplied within-model parameter priors alone do not specify these probabilities.

Let $r=(a_2,b_2)$, use the same proper prior $p(r)$ in both models, and denote their likelihoods by $L_0(r)$ and $L_1(\mu ,r)$. The unnormalized joint model targets areFor a birth move, draw an auxiliary $u$ from a positive probability density function $q(u)$ and map $(r,u)$ to $({\mu }^{\prime }=u,r^{\prime }=r)$. The reverse death move simply deletes $\mu$. This is a dimension-matching bijection with absolute Jacobian one. If the birth and death selection probabilities are $b_0$ and $d_1$, respectively, the birth and death moves for Bayesian variable selection rule givesThe reverse acceptance is the reciprocal ratio, capped at one, evaluated at the current $\mu$. Choosing $q=p_{\mu }$ cancels the added-parameter prior; a data-informed proposal may mix better but must retain its proposal factor. Within each model, also update its continuous parameters to ensure exploration.

For an explicit likelihood ratio, with $E$ and $T$ from part (ii),This supplies a complete reversible-jump Markov chain Monte Carlo move. Its acceptance probabilities enforce detailed balance between the two dimensions. The lower-dimensional state is a separate point-mass model, not a zero-probability equality test inside the continuous full-model prior.

Under the reference convention ${\alpha }_1=0$, the hypothesis becomes absence of the free group contrast ${\alpha }_2$. Introduce an inclusion indicator $J$: set ${\alpha }_2=0$ when $J=0$, and use its stated proper normal prior when $J=1$. Assign positive prior probabilities $q_0,q_1$ to those models, and retain $\mu ,{\beta }_2$ in both.

A birth draws $u$ from a proposal probability density function $g(u)$ and appends it as ${\alpha }_2^{\prime }=u$; the reverse death deletes it. The Jacobian is one. At fixed $\theta ,{\beta }_2$, the likelihood ratio isThus the birth acceptance ratio before capping iswhere $p_{\alpha }$ is the normalized $N(0,{\sigma }_1^2)$ prior. Again, taking $g=p_{\alpha }$ cancels that factor. Use the reciprocal death rule and ordinary within-model parameter updates.

After warm-up, the Markov chain Monte Carlo occupation fraction estimates the desired probability:This requires the chain to explore both models and satisfy its usual ergodic convergence conditions. Under the original continuous normal prior alone, the equality event has posterior probability zero; the explicit model indicator supplies the necessary prior atom. If intercept selection from part (v) is also performed, use two indicators and count all sampled states with $J=0$, regardless of the intercept indicator.

The parameter range is $0\le \theta \le 1$. Omitting data-only constants from the multinomial distribution likelihood, its log isFor an interior maximizer, the likelihood score must vanish:Multiply by $\theta (1-\theta )(7+4\theta )$ and collect powers to obtainThe likelihood is concave: its second derivative in the interior isTherefore an admissible interior root is the unique maximum. Endpoints must also be handled if some counts vanish. If $x_3>0$ and $x_1+x_2>0$, both endpoints have zero likelihood and the root is interior. If the maximum is at zero then $x_3=0$, and zero also satisfies the displayed polynomial; a maximum at one requires $x_1+x_2=0$, and one likewise satisfies it. Multiplying the score can introduce other endpoint roots, so not every polynomial root is an MLE: evaluate the actual likelihood on the feasible candidates.

Let $Z$ be the unobserved part of the fourth count associated with the component of probability $\theta /3$. Its complementary subcount has probability $7/12$. Conditional on the observed fourth count and a parameter value,The complete-data log-likelihood now has the simple formThe constant term can depend on the completed counts, but not on $\theta$. The split therefore converts the troublesome $\log (7+4\theta )$ term into a binomial-like maximization problem. Since $Z$ is unobserved, it should be averaged conditionally in an EM algorithm, not guessed or treated as observed.

For observed data $x$ and missing data $z$, let $p_{\theta }(x,z)$ be the complete-data model. Starting at ${\theta }^{(0)}$, the expectation-maximization algorithm alternatesThe E-step computes the conditional expectation of the complete-data log-likelihood, using the current parameter for the missing-data conditional distribution. The M-step maximizes that expected log-likelihood over a new parameter. Maximizing the likelihood of imputed mean data is equivalent only when the complete-data log-likelihood has the necessary linear dependence on the missing sufficient statistics.

For the monotonicity argument, let $q_r(z)=p_{{\theta }^{(r)}}(z\mid x)$. Jensen inequality givesTherefore an M-step increasing $Q$ cannot decrease the observed-data likelihood. Standard support and integrability conditions are understood. EM need not find a global maximum for a general model; its behavior depends on the objective and starting point.

For these counts, start with any $0<{\theta }^{(0)}<1$. The conditional missing count from part (ii) has meanThe E-step objective isIts strictly concave M-step is explicit:Equivalently, the update map is $F(\theta )=(35+360\theta )/(245+480\theta )$. Its fixed-point equation isThe unique root in $(0,1)$ is the maximum-likelihood estimateFor example, starting at $1/2$ gives $0.443299$, then approximately $0.425065$, and subsequent iterates approach the displayed root. Convergence is also transparent without relying only on the general EM theorem: $F$ is increasing andThis is positive below the feasible fixed point and negative above it. Monotonicity of $F$ prevents crossing that fixed point, so the iterates converge monotonically toward it from either side. The original log-likelihood is strictly concave and diverges to minus infinity at both endpoints, establishing that the fixed point is the unique global maximum.