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

A causal time-series representation uses only the present and past driving white noise. Thus the coefficient condition isThe series must have its stated convergence meaning. For centered white noise of positive finite variance, ${\sum }_{r\ge 0}|a_r{|}^2<\infty$ is sufficient and necessary for mean-square convergence. In the usual stable-filter convention one imposes the stronger ${\sum }_{r\ge 0}|a_r|<\infty$. A bilateral stationary linear process need not be causal: terms with $r<0$ involve future driving values.

An invertible time-series representation recovers the driving white noise from current and past observations. In the inverse series the support condition is thereforeAgain the series must converge. Stable invertibility uses ${\sum }_{r\ge 0}|b_r|<\infty$, which ensures mean-square convergence when $X_t$ has finite variance. Merely writing a bilateral inverse is not invertibility in this one-sided sense: it may require future observations. For a general correlated input $X$, square summability of $b_r$ alone is not the same sufficient condition as it is for a white noise input.

For an autoregressive moving-average model, the driving-to-output transfer function is $A(z)=\theta (z)/\varphi (z)$, and the inverse transfer function is $A(z{)}^{-1}=\varphi (z)/\theta (z)$. The causality and invertibility root criteria for an ARMA model require these respective rational functions to have power series about zero converging on a disk larger than the unit disk. If the two polynomials have no common factor, the conditions becomeIndeed, outside-disk roots leave a radius of convergence greater than one, so the coefficients decay geometrically and are absolutely summable. Conversely an uncancelled pole inside or on the unit disk prevents the required stable power series. If factors are common, apply the criterion after cancellation, to the noise-driven solution rather than additional homogeneous components. The backshift operator translates these power series into the desired one-sided filters.

The original representation has $\varphi (z)=1-2z$ and $\theta (z)=1+2z$, with roots $1/2$ and $-1/2$. There is no cancellation. Hence it is neither causal nor invertible relative to its specified driving noise. Stationarity is nevertheless possible through a two-sided solution: expanding the autoregressive inverse in negative powers givesThis is an anticausal time series representation with square-summable coefficients.

To establish the alternative representation on the same process, defineThe inverse of $1+\frac{1}{2}B$ is a stable one-sided filter. For $|z|=1$, the identities $|1+2z{|}^2=4|1+z/2{|}^2$ and $|1-2z{|}^2=4|1-z/2{|}^2$ giveThe time-series spectral density of $X$ is ${\sigma }^2|1+2e^{-i\omega }{|}^2/(2\pi |1-2e^{-i\omega }{|}^2)$; therefore the defined $\eta$ has constant time-series spectral density ${\sigma }^2/(2\pi )$. Its mean is zero, its variance is ${\sigma }^2$, and all its nonzero-lag autocovariances vanish. It is thus weak white noise. Its definition directly givesThis root reflection of an ARMA representation has roots $2,-2$, so is causal and invertible. A constant spectrum proves whiteness, not independence of non-Gaussian coordinates; no Gaussian assumption is needed for the required white noise representation.

For best linear prediction from an infinite past, let ${\mathcal{H}}_0$ be the closed linear span of $X_t$ with $t\le 0$. The causal representation expresses every such $X_t$ in present and past $\eta$ values, while invertibility puts every ${\eta }_t$ with $t\le 0$ in ${\mathcal{H}}_0$. In particular ${\eta }_1$ is orthogonal to ${\mathcal{H}}_0$. The representation at time one then gives the orthogonal projectionConsequently the best linear predictor and its error variance areOrthogonality proves optimality among linear predictors in the closed past span, without asserting that the predictor must be the conditional mean for a non-Gaussian process.

For a real weakly stationary process, extend its autocovariance by ${\gamma }_{-k}={\gamma }_k$. The exact existence of a time-series spectral density condition is that its spectral measure of a stationary time series be absolutely continuous with respect to Lebesgue measure. Absolute summability ${\gamma }_0+2{\sum }_{k\ge 1}|{\gamma }_k|<\infty$ is a useful sufficient condition, not a necessary one.

We use the conventional angular-frequency density on $[-\pi ,\pi ]$, restricted to $[0,\pi ]$ by symmetry. Under the absolute-summability condition, the Fourier series and its inverse relation areThus $2{\int }_0^{\pi }f={\gamma }_0$. If a one-sided density is instead normalized to integrate to the full variance, use $g=2f$ and omit the factor two in the inverse formula. This is the positive-frequency spectral normalization convention difference. With merely an integrable spectral density, the inverse relation remains valid; one must not assume pointwise convergence of the unweighted Fourier series. Its Fejér sums recover the density in $L^1$:

If $X$ and $Y$ are independent stationary processes, the cross covariances vanish, soLinearity of the inverse Fourier series relation, or addition of the spectral measures of a stationary time series, gives $f_Z=f_X+f_Y$. Independence can in fact be weakened to zero cross-covariances at every lag.

For the ARMA$(1,1)$ representation, away from an uncancelled unit root the transfer function givesThe usual causal case has $|\varphi |<1$; the same expression holds for the two-sided stationary solution when $|\varphi |>1$. We work in the nondegenerate case $v,w>0$ and $\varphi \ne \pm 1$; cancellations and zero-noise cases are obtained by the appropriate reduced representation or limits. Since independent white noise of variance $w$ contributes $w/(2\pi )$, putThenThe white-noise addition to an ARMA(1,1) process problem is therefore the factorizationDefineBoth are positive under the stated nondegeneracy conditions. An invertible choice isThese satisfy $|\alpha |<1$, $\lambda (1+\alpha {)}^2=P$, and $\lambda (1-\alpha {)}^2=Q$. HenceIn particular, it is the autoregressive coefficient $\varphi$ that enters the terms from the added observation noise.

To obtain a representation of the actual $Z$, defineThe stable inverse exists because $|\alpha |<1$, and spectral filtering gives $f_{\xi }=\lambda /(2\pi )$. Thus $\xi$ is weak white noise, andIts variance can equivalently be writtenWhen $\alpha \ne 0$, also $\lambda =(v\theta -w\varphi )/\alpha$; when $\alpha =0$, necessarily $C=0$ and $\lambda =A$, so the latter quotient should not be used. The model may reduce in order through cancellation. Boundary limits with $P=0$ or $Q=0$ give $|\alpha |=1$ and need not be stably invertible; the displayed stable reconstruction applies to the nondegenerate case above.

For the Box-Muller transform, draw independent $U,V$ from the uniform distribution on $(0,1)$ and setThe radial density is $r{e}^{-r^2/2}$ for $r>0$, with an independent uniform angle. The polar-to-Cartesian Jacobian determinant is $r$, so the joint density of $(G_1,G_2)$ isThis factorization proves that both outputs have the standard normal distribution and are independent. Repeating the Box-Muller transform with fresh independent uniform pairs gives independent normal outputs; discard one extra output if the desired sample size is odd.

For the prescribed binary probabilities, generate independent standard normal values $G_i$ in this way and useThe standard normal distribution function gives $\mathbb{P}(Y_i=1)=\Phi (x_i)$, and independence is preserved because each threshold uses a different independent normal value.

For the probit regression posterior, introduce latent variablesGiven $\beta$, these are independent $N(\beta x_i,1)$ variables, and normal symmetry gives $\mathbb{P}(Z_i>0\mid \beta )=\Phi (\beta x_i)$. Thus this data augmentation has exactly the observed binary likelihood. With the specified normal prior distribution, the augmented joint posterior is proportional toThe latent-normal Gibbs sampler for probit regression alternates two blocks. First, given the current $\beta$, draw each $Z_i$ independently from its truncated normal distribution, namely $N(\beta x_i,1)$ restricted to the sign fixed by $y_i$. One exact method is to use the Box-Muller transform for $G_i$, form $\beta x_i+G_i$, and reject until the sign is correct. The probability of success is positive at every finite parameter value, so the method is valid, although it can be slow for a rare sign.

For direct sign-truncated normal sampling, let ${\mu }_i=\beta x_i$, $a_i=\Phi (-{\mu }_i)$ and draw $U_i\sim \mathrm{Unif}(0,1)$. An inverse transform sampling implementation isUse suitable tail or survival-function evaluations when floating-point probabilities approach zero or one; the rejection construction remains a valid alternative.

Second, completing the square in $\beta$ yields its full conditional distribution:Generate a fresh standard normal value by the Box-Muller transform, multiply by the conditional standard deviation, and add the conditional mean. Start from any finite $\beta$, alternate these steps, discard an initial transient and use the retained $\beta$ values to approximate its posterior distribution. These are dependent Markov chain Monte Carlo samples, rather than independent posterior draws. Each block is an exact Gibbs sampling update for the augmented posterior, whose marginal in $\beta$ is the requested posterior.

The reversible-jump Markov chain Monte Carlo state comprises the model index and that model's parameter. Its unnormalized posterior density in model $j$ isFor models of equal dimension, choose model $k\ne j$ with probability $s_{jk}(\theta )$ and propose ${\theta }^{\prime }$ using density $q_{jk}({\theta }^{\prime }\mid \theta )$ on the destination parameter space. The equal-dimension reversible-jump acceptance probability isOn acceptance change both model and parameter; on rejection retain both. Choose model proposals connecting all models of positive posterior mass and within-model kernels exploring their supports, with aperiodicity, to obtain an ergodic chain. Include within-model updates targeting its conditional posterior, and use model occupation proportions after the initial transient to estimate posterior model probabilities. The formula includes the model prior probabilities, parameter prior densities, likelihoods, reverse model-selection probability, and reverse parameter-proposal density.

This direct-density formula is Metropolis–Hastings algorithm on a disjoint union of equal-dimensional spaces. If instead a move uses a deterministic bijection $T:(\theta ,u)\mapsto ({\theta }^{\prime },u^{\prime })$ with auxiliary proposal densities $g_{jk}(u\mid \theta )$ and $g_{kj}(u^{\prime }\mid {\theta }^{\prime })$, useThe Jacobian determinant is one for identity matching, but equal model dimensions alone do not make a nonlinear map volume-preserving. For a direct proposal density the change of variables is already included in that density, so no additional Jacobian factor is inserted.

For the two Poisson models, use the shape-rate convention for the gamma distribution:If $b$ denotes scale instead, replace every occurrence of the rate $b$ below by $1/b$. Introduce the inactive ${\gamma }_2$ under model $2$ with the proper pseudo-prior $\pi ({\gamma }_2)$, independent of ${\gamma }_1$. This pseudo-prior augmentation for model comparison does not change model $2$'s marginal likelihood because the inactive density integrates to one. The two augmented targets on the same positive quadrant arewhereUse a symmetric cross-model proposal that flips the model label and leaves both parameters unchanged. The matching map is the identity, with unit Jacobian determinant, and the common prior factors cancel. Hence the model-switch acceptance probabilities areOnly the second observation's factor changes because ${\gamma }_1$ is retained as its shared candidate mean. Evaluate the ratio on a logarithmic scale as ${\gamma }_2-{\gamma }_1+x_2(\log {\gamma }_1-\log {\gamma }_2)$ for the first direction.

For within-model moves, Poisson-gamma conjugacy gives exact Gibbs sampling refreshes:The two draws are conditionally independent within either model. Choose with fixed positive probabilities between this block refresh and the cross-model proposal. Both moves preserve the augmented posterior; refreshing the inactive parameter maintains its correct distribution and supports mixing. The fraction of retained states labelled $M_j$ estimates $\mathbb{P}(M_j\mid x)$.

An equivalent literal dimension-changing implementation deletes ${\gamma }_2$ on a proposed $1\to 2$ move, and on a $2\to 1$ move draws an auxiliary $u\sim \mathrm{Gamma}(a,b)$ and sets ${\gamma }_2=u$. The birth matching has $1+1=2+0$ dimensions and unit Jacobian determinant. Its proposal density cancels the new parameter's prior density, giving the same acceptance ratios above for symmetric selection of move directions. This provides suitable moves without storing an inactive parameter.

For an independent check of Poisson mean equality model comparison, the exact Bayes factor isIt follows by integrating the two gamma likelihood kernels; the common factorial terms cancel. With equal model priors, the exact posterior probability of $M_2$ is $m_2/(m_1+m_2)$, which can also be used to check the RJ-MCMC occupation estimate.