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

For a series with constant mean $\mu$, write $X_t=x_t-\mu$. An autoregressive moving-average process has the formwhere $w_t$ is zero-mean white noise with variance ${\sigma }_w^2>0$. The backshift operator is $B{X}_t=X_{t-1}$. Define the autoregressive polynomial and the moving-average polynomial byThe autoregressive operator and moving-average operator give the concise equationOrders are normally specified minimally, with nonzero last coefficients and no common polynomial factor. A displayed unreduced equation can have larger orders than the process actually needs.

Causality means that the stationary process can be recovered from present and past noise:Invertibility means that the noise can be recovered from present and past observations:The causality and invertibility root criteria for an ARMA model, for the reduced representation, areIndeed the transfer series are $\sum {\psi }_j{z}^j=\Theta (z)/\Phi (z)$ and $\sum {\pi }_j{z}^j=\Phi (z)/\Theta (z)$. When all denominator roots lie outside the closed unit disc, each rational function is analytic on a larger disc, so its Taylor coefficients decay geometrically and are absolutely summable. Conversely, a genuine denominator zero in the closed unit disc gives a pole incompatible with such a convergent filter. The reduced-representation proviso prevents a cancelled root from being mistaken for an obstruction.

For order identification by autocorrelation cutoffs, first estimate the sample autocorrelation function and sample partial autocorrelation function from a suitably stationary, centred series. The autocorrelation function is $\rho (h)=\mathrm{Cov}(X_t,X_{t-h})/\mathrm{Var}(X_t)$. The lag-$h$ partial autocorrelation function is the last coefficient in the best linear regression of $X_t$ on $X_{t-1},\dots ,X_{t-h}$, equivalently the correlation left after removing the intervening lags.

For a pure moving-average model of order $q$, observations more than $q$ lags apart share no noise, so the population ACF is zero beyond $q$, while the PACF usually tails off. For a pure autoregressive model of order $p$, the PACF is zero beyond $p$, while the ACF usually decays, possibly with damped oscillation. Thus a clear sample ACF cutoff suggests an MA order and a clear PACF cutoff suggests an AR order. For a mixed ARMA process both usually tail off, so these plots guide a small set of candidate orders rather than determine $p,q$ uniquely. Sampling variation blurs cutoffs; compare fitted candidates and check that their residual ACF resembles white noise. Common-factor cancellation can also disguise the orders in a displayed equation.

The displayed equation has nominal orders $p=1,q=2$, withThe AR root is $2$, and the MA roots are $2$ and $10/9$, all strictly outside the unit disc. Thus its stationary solution is causal and invertible. There is also common factor cancellation in an ARMA model: the convergent inverse of $1-0.5B$ cancels the shared factor, givingConsequently the displayed representation is ARMA(1,2), but its minimal order is ARMA(0,1). This also follows without formal operator cancellation: the difference between any solution and the right-hand side satisfies $d_t=0.5d_{t-1}$. A stationary finite-variance difference must have zero variance, since stationarity would give $\mathrm{Var}(d_t)=0.25\mathrm{Var}(d_t)$; its constant mean must also vanish. A transient chosen by an arbitrary initial condition is not the stationary process whose ACF is requested.

For this moving-average model, independence of the noise giveswhile $\gamma (h)=0$ for $|h|\ge 2$, since the two noise sets are disjoint. Therefore its autocorrelation function isThe reduced moving-average inverse is $w_t={\sum }_{j\ge 0}0.9^j{x}_{t-j}$, whose coefficients are absolutely summable, confirming invertibility directly.

Use the standard linear Gaussian state-space model convention: the initial state is independent of both mutually independent iid noise sequences. This independence from the initial state is needed for the usual forecast recursion. Let ${\mathcal{Y}}_t=\sigma (y_1,\dots ,y_t)$ and write $m_{t-1}=x_{t-1}^{t-1}$ and $C_{t-1}=P_{t-1}^{t-1}$. Taking conditional expectations in the state equation givesTo justify the covariance, the forecast error is $\Phi (x_{t-1}-m_{t-1})+w_t$. Its two terms have zero cross-covariance because $w_t$ is independent of the previous state and observations. Their covariances therefore add. In this Gaussian model the conditional error covariance depends on the parameters and time but not on the realized data, so it also equals the unconditional forecast-error covariance used in the notation of the question.

For completeness, derive the filtering step rather than leave the algorithm unspecified. Set $m_t^{-}=x_t^{t-1}$ and $C_t^{-}=P_t^{t-1}$. Conditionally on the previous data, the observation has mean $A{m}_t^{-}$ and covarianceDefine the innovation process and gain byThe state forecast error and innovation are jointly normal, with cross-covariance $C_t^{-}{A}^{\top }$. The residual $x_t-m_t^{-}-K_t{ϵ}_t$ is uncorrelated with ${ϵ}_t$ by the definition of $K_t$, hence independent of it by joint normality. Its mean is zero and its covariance is $C_t^{-}-C_t^{-}{A}^{\top }{\Sigma }_t^{-1}A{C}_t^{-}$. This proves the conditional multivariate normal distribution updateInitialize $x_0^0={\mu }_0,P_0^0={\Sigma }_0$. For $t=1,\dots ,n$, predict using the first boxed pair, store that forecast, form ${ϵ}_t,{\Sigma }_t,K_t$, and update using the second boxed pair. The updated state then starts the next iteration. This is the Kalman filter, and supplies the entire forecast sequence in one forward pass.

The requested unconditional innovation moments areIn fact ${ϵ}_t\mid {\mathcal{Y}}_{t-1}\sim N_q(0,{\Sigma }_t)$, with ${\Sigma }_t$ deterministic for fixed parameters. Thus the innovation is independent of the past observations, and hence the innovations are mutually independent. Equivalently, orthogonality to the past together with joint normality gives this independence.

The chain rule for densities now yields the Gaussian innovation likelihood, treating the specified initial mean and covariance as known:Both ${ϵ}_t$ and ${\Sigma }_t$ are functions of the trial parameters, so the filtering recursion must be run again for each candidate $\Theta$. The ordinary density formula assumes ${\Sigma }_t$ positive definite; for example positive definite $R$ ensures this. Singular observation laws instead require densities on their Gaussian support.

Evaluate the log-likelihood numerically by the forward recursion, and maximize it over admissible $\Phi ,A,Q,R$, enforcing the covariance constraints. A Cholesky decomposition ${\Sigma }_t=L_t{L}_t^{\top }$ gives $\log |{\Sigma }_t|=2{\sum }_j\log (L_t{)}_{jj}$ and the quadratic term $\Vert L_t^{-1}{ϵ}_t{\Vert }^2$, avoiding an explicit matrix inverse. Covariance parameterizations or constrained optimization keep $Q,R$ positive semidefinite and the innovation covariance nonsingular. An attained maximizer is a maximum-likelihood estimator; in general this is a numerical optimization problem and neither uniqueness nor global convergence of a local optimizer is automatic.

Choose a proposal probability density function $g$ and a finite envelope constant $M$ such that $f(y)\le Mg(y)$ wherever the target has positive density. A rejection sampling trial draws $Y\sim g$ and, independently, $U\sim \mathrm{Uniform}(0,1)$; accept $Y$ whenotherwise repeat with fresh independent draws. For any measurable set $D$,In particular $\mathbb{P}(\mathrm{accept})=1/M$, so the conditional accepted density is $f$. More explicitly, summing over the possible first accepted trial givesThus the algorithm terminates almost surely and its output has exactly the requested law.

The envelope constant satisfies $M\ge 1$ and is the expected number of proposal trials for one output, by the geometric distribution of the waiting time. This is the quantity often called the rejection rate $M$. The actual probability of rejecting a particular proposal is $1-1/M$, and the expected number of rejections before one acceptance is $M-1$; these conventions should not be conflated.

To generate a standard normal distribution, first generate its absolute value, whose half-normal distribution density isWith the given exponential distribution proposal, the density ratio isIts maximum is at $y=\lambda$, hence the normal rejection sampling with an exponential envelope constant isDraw $Y\sim \mathrm{Exp}(\lambda )$ and a fresh uniform $U$; accept when $U\le e^{-(Y-\lambda {)}^2/2}$. An accepted $Y$ has density $h$. Independently choose a sign, positive or negative with probability $1/2$, and output $Z$ with that sign and magnitude $Y$. On either half-line its density is $h(|z|)/2=(2\pi {)}^{-1/2}{e}^{-z^2/2}$, proving the normal output law.

To optimize the rejection rate, differentiateIts derivative is $\lambda -1/\lambda$, negative before one and positive after one; equivalently its second derivative is $1+1/{\lambda }^2>0$. ThusThe rejection probability is consequently about $0.23983$.

To use only uniforms, apply inverse transform sampling: for an independent $U_1\in (0,1)$,has survival probability $\mathbb{P}(Y>y)=e^{-\lambda y}$ for $y\ge 0$. Use a second independent uniform for acceptance, and, after acceptance, a third independent uniform to choose the sign. Repeat with fresh uniforms after any rejection. No exponential generator is then needed.

An unordered bootstrap sample is determined by the multiplicities $N_i\ge 0$ of the $n$ distinct observations, with ${\sum }_i{N}_i=n$. Conversely each such vector specifies exactly one sample up to rearrangement. Place $n$ identical marks into $n$ boxes using $n-1$ separators. The stars and bars argument therefore givesThese bootstrap count vectors are not equally likely: the probability of a vector is $n!/(n^n{\prod }_i{N}_i!)$, counting its possible ordered sequences. The ordered index sequences themselves number $n^n$.

The first simulation uses independent standard Cauchy distributions and returnsThe tail probability and estimator moments areThe indicator count has a binomial distribution with parameters $n,\theta$; the infinite variance of a Cauchy observation is irrelevant because this estimator averages bounded indicators.

For the next code block, each row is an independent sample of size $n$ from the empirical distribution of the original $x_i$. Its entry in the vector $v$ is the corresponding bootstrap indicator mean ${\hat{\theta }}_b^{\ast }$. The last expression is the variance of these $B=199$ replicate estimates, with denominator $B-1$:Condition on the original data and write a star on bootstrap expectations. Every resampled indicator has a Bernoulli distribution with success probability $\hat{\theta }$, and the $n$ draws within a row are independent. Thus the conditional bootstrap variance of a sample mean givesIt estimates the sampling variance of the tail-probability estimator, not the variance of the raw Cauchy data. Its square root is an estimated standard error. Finite $B$ leaves bootstrap simulation noise. Averaging over the original data, its expectation is $(n-1)\theta (1-\theta )/n^2$, so it is not exactly unbiased for the true variance $\theta (1-\theta )/n$ at finite $n$.

In the importance-sampling block, the transformation $Y=2/(1-U)$ gives, for $y\ge 2$,The variable $Y$ is a proposal draw supported entirely in the integration region. The deterministic value computed from it is the importance sampling weightConsequentlyThe code uses ordinary importance sampling of a Cauchy tail, not self-normalization; $Y$ itself is neither an observation from the target Cauchy law nor the estimator's summand. The weight compensates for sampling from its proposal law.

The final block bootstraps these same proposal observations $Y_i$ and recomputes the weight average in each row. Since $W$ is a deterministic function of $Y$, this is equivalent to resampling the observed weights $W_i$. If ${\tilde{\theta }}_b^{\ast }$ denotes a replicate, the final output is its sample variance over the $B$ replicates. Its conditional expectation isIts expectation over the original proposal sample is $(n-1)\mathrm{Var}(W)/n^2$.

The variance reduction can be calculated, rather than inferred just from the appearance of the weights. Write $A=\arctan (1/2)$. Thenbecause an antiderivative of the integrand without its prefactor is $\frac{1}{2}(\arctan y-y/(1+y^2))$. HenceBy comparison, $\theta (1-\theta )\approx 0.125803$. For $n=100$ the true estimator variances are approximately $9.55253\times 10^{-7}$ and $1.25803\times 10^{-3}$ respectively. The importance estimator has about 1/1,317 of the indicator estimator's variance, because all proposals land in the tail and the weights remain in the narrow interval $[2/(5\pi ),1/(2\pi ))$.

Thus the last bootstrap output should usually be much smaller than the first, and its expected value is smaller by the same variance ratio. This is not a deterministic ordering of two finite-run outputs: if all the original Cauchy draws happen to be at most two, the first bootstrap variance is zero, while unequal importance weights can still give a positive second variance.

The Metropolis–Hastings algorithm proposes $y\sim q(\cdot \mid x)$ and moves from the current state $x$ with acceptance probabilityOn rejection it retains $x$. The target need only be known up to its normalizing constant, but it must define a proper probability distribution. The accepted off-diagonal probability flow iswhich is symmetric in $x,y$ and proves detailed balance. This permits considerable flexibility in proposals, including joint moves when a posterior distribution has strongly correlated parameters, but requires proposal tuning and can waste computation on rejections.

A Gibbs sampler instead updates a coordinate or block from its exact full conditional distribution, holding the other coordinates fixed. This is a special Metropolis–Hastings algorithm update with acceptance probability one: the joint density factors into the unchanged marginal density and the proposed conditional density, which cancel in the acceptance ratio. The coordinate update preserves the target; composing such updates therefore preserves it too, although an entire deterministic sweep need not itself be reversible.

Use Gibbs sampling when the conditional laws are easy to simulate and suitably chosen blocks mix well. Scalar updates can mix poorly under strong posterior distribution correlation, whereas a blocked Gibbs sampler can alleviate that dependence. Prefer a flexible Metropolis proposal when conditional draws are difficult or a well-tuned joint proposal is more effective. Variations include random or systematic coordinate scans, random-walk Metropolis algorithms, independence proposals, and using Metropolis updates for difficult blocks inside a Gibbs scheme. Both methods generate dependent draws; usable long-run estimates require an appropriate irreducible, aperiodic chain and assessment of mixing, not merely a formal update rule.

For the change-point calculation, put $v={\sigma }^2$ andMultiplication of the Gaussian likelihood function by the specified prior distribution gives the formal posterior distribution kernel, with respect to $d\theta d\varphi dv$ and counting measure on $m$,

There is an empty component under an improper prior at $m=n$. The sum $S_n$ is independent of $\varphi$, whose prior distribution is flat on all of $\mathbb{R}$. For every fixed finite $\theta$ and $v>0$, the kernel is a strictly positive constant in $\varphi$, henceThe configuration has positive prior probability $1/n$. Integrating over, for example, a bounded interval of $\theta$ and $1\le v\le 2$ already makes the joint normalization infinite. Thus the printed prior distribution does not produce a posterior probability distribution. This failure is present for every data set, not just for a special arrangement of observations.

The formal full conditionals expose exactly where the requested sampler breaks. For $m<n$, completing the two mean squares givesThe variance kernel is that of $\mathrm{IG}(n/2,S_m/2)$ when $S_m>0$, and the discrete conditional probabilities are proportional to $e^{-S_m/(2v)}$. At $m=n$ the conditional kernel for $\varphi$ is constant and cannot be normalized; it is not a normal distribution with an infinite variance. Therefore no valid Gibbs sampler, and no value of $P(m=n\mid x)$, exists for the model as printed. A chain that simply substitutes an arbitrary distribution for this missing conditional would be sampling a different model.

Here is a fully specified proper-prior repair retaining the no-change configuration. Choose fixed hyperparametersindependently, with $t_{\theta },t_{\varphi },a_0,b_0>0$. Its Gaussian change-point posterior with proper priors isIt is normalizable even for degenerate data: the likelihood function is at most $(2\pi {)}^{-n/2}{v}^{-n/2}$, and its integral against the variance prior distribution is bounded byThe normal mean priors integrate to one, and the normalizing constant is positive. This proves posterior propriety for the repaired model.

Completing the conditional squares gives the Gaussian change-point Gibbs updates. DefineThe four updates areAt $m=n$, $M_{\varphi }={\mu }_{\varphi },V_{\varphi }=t_{\varphi }$, giving its proper prior distribution unchanged. Starting from any finite means, positive variance and allowed $m$, cycle through these four updates using the latest values. The sum-of-squares terms can be evaluated from cumulative data sums; normalize the discrete log weights after subtracting their maximum to avoid numerical underflow. The proper positive conditionals allow movement among all split values and throughout the continuous parameter support.

Under this specified repair, $P(m=n\mid x)$ is the posterior probability of no change within the observed sequence: all $n$ data points belong to the first mean regime. For $N$ retained draws after initialization has ceased to affect the estimates, its Markov chain Monte Carlo estimator isAlternatively average the displayed conditional probability of $m=n$ at the sampled continuous parameters, giving an estimate by Rao-Blackwellization. Autocorrelation must be accounted for in its Monte Carlo uncertainty. These are probabilities and estimators for the explicitly repaired model; the original unnormalizable kernel supplies neither.

For the expectation-maximization algorithm, let $L(\vartheta )=\int f(x,z;\vartheta )dz$. At iteration $t$, the E-step forms the conditional distribution of the missing data at the current parameter and the functionThe M-step isIn this maximization the conditional distribution defining the expectation is held at the old parameter; it is not recomputed as $\vartheta$ varies. The Jensen inequality givesby writing the likelihood function ratio as a conditional expectation of $f(x,Z;\vartheta )/f(x,Z;{\vartheta }^{(t)})$. Hence an exact M-step does not decrease the observed-data likelihood function. The algorithm does not in general guarantee convergence to its global maximum; initialization and the likelihood function's geometry matter.

For the specified diagonal bivariate normal distribution, use $v_j={\sigma }_j^2>0$. The two coordinates are independent, as are the four observations. Therefore the known second coordinate of the fourth observation gives no information about its first coordinate conditional on the parameters. The E-step isThese two moments suffice for the complete Gaussian log-likelihood. Discarding terms independent of the proposed parameters, its conditional expectation isDifferentiating with respect to each mean gives four times that mean minus the corresponding expected sum of observations. Thus the updates in EM for an independent missing normal coordinate areFor either variance the remaining criterion is $-2\log v-S/(2v)$; its derivative vanishes at $v=S/4$, which is its maximum when $S>0$. Evaluate the expected residual sum at the new means to obtainThe missing contribution isIts first term is essential. Substituting only the imputed mean for $Z$ would omit its uncertainty and would not maximize the required expected complete log-likelihood.

As a check, the observed-data likelihood function integrates out the missing first coordinate, leaving three first-coordinate and four second-coordinate normal densities. Put ${\overline{x}}_1={\sum }_{i=1}^3{x}_{i1}/3$ and $s_1^2={\sum }_{i=1}^3(x_{i1}-{\overline{x}}_1{)}^2/3$. The first-coordinate recurrence impliesThus the iterates converge to the three-observation mean and maximum-likelihood variance, while the second-coordinate estimates are reached after one update. If the observed values in a coordinate are all equal, its likelihood function instead has a zero-variance boundary degeneracy, and there is no positive-variance maximizer; the variance formulas then describe that boundary limit. For nonzero empirical spreads the displayed updates give the ordinary nonsingular EM iteration.