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

Series 1 wanders over a changing level rather than fluctuating around a stable local mean. Its sample autocorrelation function is strongly positive and decreases very slowly. This is the usual diagnostic evidence for an ordinary unit root: an autoregressive polynomial containing $1-z$, with a zero at $z=1$, and a stationary model after first differencing. The plots support an integrated model, rather than specifying the number of its remaining stationary autoregressive or moving-average terms.

Series 2 has a pronounced oscillation with period about six observations. Its sample autocorrelation alternates between large positive and negative values with little damping: approximately positive at multiples of six and negative halfway between. Together with the changing amplitude, this suggests a conjugate pair of unit-circle zeros nearThe associated real autoregressive factor is $1-2\cos (\pi /3)z+z^2=1-z+z^2$. A targeted filter $1-B+B^2$ removes this pair; the broader seasonal difference operator $1-B^6$ also contains it but introduces additional differencing factors. This is the oscillatory unit-root diagnosis from an undamped sample autocorrelation.

Thus Series 1 suggests a zero at 1; Series 2 suggests a conjugate pair on the unit circle at a seasonal frequency. These are model diagnoses, not deductions of exact roots from a finite sample. A stationary model very close to a unit root can look similar, and an undamped periodic covariance can also arise from a stationary random sinusoid. The figure does not identify exact orders or prove nonstationarity by itself.

Among the supplied candidates, choose ARMA(2,1). It has the smallest Akaike information criterion. Relative to ARMA(2,2), its AIC improvement is $1.71$; the additional second moving-average estimate is only half a standard error from zero and increases the log likelihood function by only about $0.14$. Dropping it is a reasonable parsimony choice, although that AIC gap is small. ARMA(1,1) has AIC larger by $17.55$, a much clearer loss of fit. In the chosen model the second autoregressive term is about seven standard errors from zero, so it should not be dropped merely to obtain order one. The first autoregressive term is less precisely estimated; this does not justify automatically deleting it without fitting and comparing the reduced candidate.

Using the usual positive-sign moving-average convention, the fitted autoregressive moving-average model isThese are plug-in values rounded as given, not exact population parameters. The model has zero mean. Its polynomials are $\varphi (z)=1+0.30z-0.28z^2=(1-0.4z)(1+0.7z)$ and $\theta (z)=1+0.5z$.

With angular frequency $\omega \in [-\pi ,\pi ]$, so that the autocovariance is $\gamma (h)={\int }_{-\pi }^{\pi }{e}^{ih\omega }{f}_X(\omega )d\omega$, the spectral density of a stationary process isThe frequency convention makes the normalization unambiguous.

The autoregressive polynomial factors as $(1-0.4z)(1+0.7z)$, with zerosBoth have modulus greater than one. The causality root criterion for an autoregressive model therefore gives a causal stationary solution. The moving-average polynomial has its only zero at $-2$, also outside the unit circle, so the invertibility of a moving-average model holds. There is no common root to cancel.

To use unit-variance white noise, put $W_t=Z_t/2$. One suitable pair isThen $\tilde{\varphi }(B)X_t=\tilde{\theta }(B)W_t$ with $W\sim \mathrm{WN}(0,1)$. The factor two changes the innovation scale, not the zero of the moving-average polynomial.

Write $X_t={\sum }_{j\ge 0}{\psi }_j{Z}_{t-j}$, with innovation variance four. The generating function of this linear process isEquating coefficients gives ${\psi }_0=1$, ${\psi }_1=0.2$ and ${\psi }_j=-0.3{\psi }_{j-1}+0.28{\psi }_{j-2}$ for $j\ge 2$. Thus the first five coefficients, including lag zero, areFor all remaining lags, partial fractions giveThe coefficients are absolutely summable, so the series converges in mean square and defines the causal moving-average expansion. In the unit-variance convention of part (c), $X_t={\sum }_{j\ge 0}{c}_j{W}_{t-j}$ with $c_j=2{\psi }_j$; its first five coefficients are $(2,0.4,0.44,-0.02,0.1292)$. This is the two-geometric-coefficient expansion of a causal ARMA(2,1) process.

White noise orthogonality gives, for any integer $h$,To see this directly, expand the covariance of the two convergent series. Only matching noise indices contribute. Cauchy-Schwarz inequality makes the coefficient-product sum finite, justifying the covariance limit.

There is also a closed expression. Set $a=9/11$, $d=2/11$, $r=0.4$, $s=-0.7$. For $h\ge 0$,and negative lags follow by symmetry. Dividing by the expression at zero gives the same autocorrelation function.

Interpret stationarity in the usual second-order time-series sense and assume nondegenerate noise, ${\sigma }_z^2>0$. For a two-sided autoregressive equation, the missing existence condition isIt is important to separate this from causality. If $|\varphi |<1$, the unique stationary solution is $X_t={\sum }_{j\ge 0}{\varphi }^j{Z}_{t-j}$. If $|\varphi |>1$, there is still a stationary solution, but it is anticausal:Both expansions converge in L2 because their coefficients are square summable. Substitution verifies the equation. Their means are zero and their covariance functions depend only on lag. Uniqueness follows by iterating the equation backward in the first case and forward in the second: the remainders ${\varphi }^n{X}_{t-n}$ or ${\varphi }^{-n}{X}_{t+n}$ tend to zero in L2 for any stationary finite-variance solution. This is the stationary versus causal solution of a two-sided AR(1) equation.

For $\varphi =\pm 1$, iteration givesThe variance of the right side is $n{\sigma }_z^2$. The variance of the left side is at most $4\mathrm{Var}(X_t)$ by stationarity and Cauchy-Schwarz inequality. These are incompatible as $n\to \infty$. Thus no weakly stationary finite-variance solution exists at those unit roots.

If the intended claim includes a causal innovation representation, its condition is instead $|\varphi |<1$, as in the next part. The stated white noise equation alone does not say that $Z_t$ is orthogonal to the past of $X$. If zero innovation variance is allowed, the unit-root exclusion has degenerate exceptions, such as random constant solutions when $\varphi =1$; the nondegenerate convention is necessary for the asserted nonexistence.

The condition for a causal linear-filter solution is $|\varphi |<1$. Its mean-square expansion is $X_t={\sum }_{j\ge 0}{\varphi }^j{Z}_{t-j}$. Summing the matching white noise terms givesThe assumed orthogonality of every $W_s$ to every $Z_t$ extends to every $X_t$ by L2 convergence of that expansion. Hence the added-noise process has mean zero andThis depends only on lag, proving weak stationarity. Its autocorrelation has the same geometric tail as the latent autoregression, but its positive-lag correlations are reduced by the additional variance at lag zero. This is the autocovariance of an AR(1) process observed with white noise. Strict stationarity or Gaussianity is not implied by white noise covariance assumptions alone.

Apply $1-\varphi B$ to the observed process and call the result $V_t$. ThenIts only nonzero covariance lags areSeek an invertible moving-average factor $V_t={\epsilon }_t+\vartheta {\epsilon }_{t-1}$ with $\mathrm{Var}({\epsilon }_t)=\nu$. Matching these covariances requires $\nu (1+{\vartheta }^2)=A$ and $\nu \vartheta =C$. Solving givesThese are the three requested parameters in the positive-sign moving-average convention. For nonzero total noise variance,so $|\vartheta |<1$ and the larger quadratic root gives the invertible factor.

Covariance matching alone would not identify arbitrary processes in distribution. To obtain an actual representation on the given space, defineThe series converges in L2. The spectrum of $V$ is $(A+2C\cos \omega )/(2\pi )=\nu |1+\vartheta e^{-i\omega }{|}^2/(2\pi )$, so this filtered process has constant spectrum $\nu /(2\pi )$ and is white noise. ThusThis is the invertible ARMA factorization of an AR(1)-plus-noise process. If ${\sigma }_w^2=0$, it reduces to AR(1); if $\varphi =0$, it reduces to white noise. If ${\sigma }_z^2=0$ and ${\sigma }_w^2>0$, then $\vartheta =-\varphi$, the common factor cancels and $Y=W$. Therefore the orders are at most (1,1); no unnecessarily minimal-order claim is made in those degenerate cases.

Take the latent autoregression as the scalar state $S_t$. A state-space model isThe transition and observation matrices are both scalar, $F=\varphi$, $H=1$. State-noise variance is $Q={\sigma }_z^2$, observation-noise variance is $R={\sigma }_w^2$, and the cross-noise covariance is zero at every pair of times.

A complete stationary initialization isIt is orthogonal to future state noise and to all observation noise. This is the stationary initialization of a scalar linear state-space model. This specifies the initial state in terms of the actual given two-sided white noise sequence, as well as its second-order law; simply starting from zero would give transient rather than stationary observations.

If a Gaussian state-space specification is intended, the complete specialization is $S_0\sim N(0,Q/(1-{\varphi }^2))$, independent of the future iid Gaussian state and observation noises, themselves independent with variances $Q,R$. Under the printed assumptions alone, Gaussian distributions and independence cannot be deduced from white noise orthogonality; the equations and stationary-series initialization above give the exact second-order representation without adding them.

Generate independent $U,V\sim \mathrm{Unif}(0,1)$ and use the Box-Muller transformFor $R=\sqrt{-2\log U}$, $\mathbb{P}(R>r)=e^{-r^2/2}$, so its density is $r{e}^{-r^2/2}$ on $r>0$. The angle $2\pi V$ is uniform on $[0,2\pi )$ and independent of $R$. Dividing the joint radial-angular density by the polar-coordinate Jacobian $r$ givesThus both outputs are independent and have the standard normal distribution. Independent pairs of uniforms give further independent normal observations. The null event $U=0$ is excluded in implementation so the logarithm is finite.

For fixed parameters $(a,b)$, every new error is independent of the previous observations. Put ${\phi }_0(z)=(2\pi {)}^{-1/2}{e}^{-z^2/2}$. Since $X_0=X_1=0$, the first three conditional densities areandIn general,Its density is ${\phi }_0(x_{t+2}-a{x}_{t+1}-b{x}_t)$. The recursion is interpreted from $t=0$, as required to define $X_2$ from the two given initial values. These are parameter-conditional sampling distributions; integrating over the prior would instead give predictive mixtures.

Use the chain rule for conditional densities. With $x_0=x_1=0$, the likelihood function of the nondegenerate observations isThe first residual is $x_2$ and carries no parameter information. There is no ordinary joint Lebesgue density for $X_1,\dots ,X_n$ because $X_1=0$ deterministically; this expression is the conditional likelihood function given the fixed initial values, or the likelihood function on $x_2,\dots ,x_n$. The initial point mass is parameter independent and has no effect on inference. This is the conditional likelihood of an initialized Gaussian AR(2) process.

Define the sufficient data sums, all over $0\le t\le n-2$,Multiplying the likelihood function by the independent standard-normal priors and collecting the quadratic terms givesLet $\Lambda =\left(\begin{array}{cc}A& C\\ C& D\end{array}\right)$ and $\Delta =AD-C^2$. This precision matrix is $I+\sum v_t{v}_t^T$, where $v_t=(x_{t+1},x_t{)}^T$, so it is positive definite even for a short or singular design. Completing the square provesThe fully normalized posterior density isThis is Gaussian conjugacy for an initialized AR(2) regression.

Completing each one-dimensional square, or using the supplied conditional-normal identity, givesWhen $n=2$, all regressor sums vanish and the posterior remains the independent standard-normal prior. No stationary-parameter restriction is imposed: the specified prior is on all of ${\mathbb{R}}^2$, and the finite initialized chain is defined for all $(a,b)$.

Start from any $(a^{(0)},b^{(0)})$. For independent draws $G_{1,k},G_{2,k}$ from the standard normal distribution at each sweep, the Gibbs sampler can be implemented asThe second update must use the newly drawn first coordinate. The transform in part (a) supplies the needed independent Gaussian draws. Each conditional update preserves the joint posterior, so their composition does too.

Convergence is particularly transparent here. Centering at the posterior mean givesPositive definiteness gives $C^2/(AD)<1$, so this is a stable Gaussian autoregression. The full sampler has the desired posterior as its limiting invariant law. This is the linear contraction of a two-coordinate Gaussian Gibbs sweep.

For a posterior-integrable function $\varphi$, its posterior expectation is estimated after a burn-in $K$ byThe ergodic theorem justifies this Monte Carlo average. The retained values of $\varphi$ also approximate its posterior distribution and quantiles. If a Bayesian point estimate is requested under squared-error loss, the posterior mean is the appropriate estimate, when its needed moments exist. An arbitrary function need not have a finite posterior mean; integrability must be assumed for the displayed target. Successive Gibbs draws are correlated, so uncertainty in the Monte Carlo average should use chain-aware error estimates rather than treating them as iid.

Let $w(x)=f(x)/g(x)$ where $g(x)>0$, and set it to zero on the g-null set where $g=0$. The domination assumption makes $f=0$ there and gives $0\le w\le M$. Integrating it also gives $M\ge 1$. For iid proposals $X_i\sim g$, the importance sampling estimator isCauchy-Schwarz inequality under $f$ gives $|\theta |\le (\int {\varphi }^{2}f{)}^{1/2}<\infty$. Direct integration proves unbiasedness, ${\mathbb{E}}_g[\varphi (X)w(X)]=\theta$. MoreoverThusThe iid central limit theorem applies to these finite-variance weighted observations:If the variance is zero, this denotes the point mass at zero. This is the bounded-weight importance-sampling moment bound.

For each proposal $X\sim g$, independently draw $U\sim \mathrm{Unif}(0,1)$ and accept it whenThe ratio is at most one, as required. For any measurable set $A$,so the acceptance probability is $p=1/M$ and the conditional distribution of an accepted proposal has density $f$. Repeating independent trials until acceptance therefore gives an exact draw from $f$, and repeating the whole procedure gives iid target draws. This proves rejection sampling.

The number of proposals for one successful draw is geometric with mean $M$. Thus $M$ is also the expected proposal cost of this exact simulation method.

Apply independent uniforms to the $n$ proposals, and let $I_i$ indicate acceptance. Each $I_i$ is Bernoulli with success probability $p=1/M$, and the indicators are independent because the pairs $(X_i,U_i)$ are independent. HenceFor a fixed acceptance pattern, the values at accepted positions are independent with density $f$, by the single-trial conditional density calculation. The same product law holds for every pattern of a given size. Consequently, conditional on $N=m$, the ordered accepted observations have joint density ${\prod }_{j=1}^{m}f(y_j)$. This is the binomial count and iid values in fixed-budget rejection sampling; the count provides no information about those target values.

There is a finite-sample empty-output event, with probability $(1-1/M{)}^n$. Any estimator dividing by $N$ must be defined separately on $N=0$.

Write $p=1/M$ and $v_f=\int {\varphi }^{2}f-{\theta }^2$. The Bernoulli central limit theorem givesIf $M=1$, all proposals are accepted and this limit is degenerate.

Conditional on $N=m$, the accepted sample is iid from $f$. Since $N/n\to p>0$, its size tends to infinity, and the ordinary target-sample CLT therefore givesAssign any fixed value to the estimator on $N=0$; the probability of that event tends to zero, so it does not change this limit. This is the random-count central limit theorem for accepted rejection samples.

For comparison at the same budget of $n$ proposals, Slutsky's theorem rescales this asThe distinction between accepted-sample size and proposal count is essential for the final efficiency comparison.

At a common proposal budget, the two asymptotic variances areThus prefer the estimator with the smaller proposal-budget variance; the stated assumptions do not give a universal winner. Importance sampling uses every proposal and avoids an empty accepted sample, but those facts alone do not establish variance dominance. The bound from part (a) only says $v_{\mathrm{IS}}\le M{v}_f+(M-1){\theta }^2$. If $\theta =0$, it does imply that importance sampling is at least as efficient asymptotically.

For explicit nonconstant counterexamples in both directions, take $g(x)=1$ and $f(x)=2x$ on $0<x<1$, and zero outside, with the sharp envelope $M=2$. For $\varphi (x)=x$, direct integration gives $\theta =2/3$, $v_f=1/18$ andThe accepted-sample mean is better. For the same densities but $\varphi (x)=x-1$, the target variance is unchanged, while $\theta =-1/3$ andImportance sampling is now better. This is the variance reversal under additive shifts of an importance-sampling integrand.

There is a useful distinction about the Rao-Blackwell theorem. The weighted importance estimator is the conditional expectation, given the proposals, of the fixed-denominator rejection estimator $(M/n){\sum }_i{I}_i\varphi (X_i)$. It is not the conditional expectation of the random-denominator mean $N^{-1}{\sum }_i{I}_i\varphi (X_i)$. Rao-Blackwell variance reduction for the former therefore does not prove a comparison with the latter. This is the Rao-Blackwell identity for a fixed-denominator rejection estimator. The appropriate answer is the proposal-budget variance comparison of importance and rejection sampling, with the empty-output convention specified.