The SIM2D Procedure

Theoretical Development

It is a simple matter to produce an $N(0,1)$ random number, and by stacking k $N(0,1)$ random numbers in a column vector, you can obtain a vector with independent standard normal components $\mb {W} \sim N_ k(\Strong{0},\mb {I})$. The meaning of the terms independence and randomness in the context of a deterministic algorithm required for the generation of these numbers is subtle; see Knuth (1981, Chapter 3) for details.

Rather than $\mb {W} \sim N_ k(\Strong{0},\mb {I})$, what is required is the generation of a vector $\mb {Z} \sim N_ k(\Strong{0}, \mb {C})$—that is,

\[  \mb {Z} = \left[ \begin{array}{c} Z_1\\ Z_2\\ \vdots \\ Z_ k \end{array} \right]  \]

with covariance matrix

\[  \mb {C} = \left( \begin{array}{cccc} C_{11} &  C_{12} &  \cdots &  C_{1k}\\ C_{21} &  C_{22} &  \cdots &  C_{2k}\\ &  \ddots & & \\ C_{k1} &  C_{k2} &  \cdots &  C_{kk} \end{array} \right)  \]

If the covariance matrix is symmetric and positive definite, it has a Cholesky root $\mb {L}$ such that $\mb {C}$ can be factored as

\[  \mb {C}=\mb {LL}’  \]

where $\mb {L}$ is lower triangular. See Ralston and Rabinowitz (1978, Chapter 9, Section 3-3) for details. This vector $\mb {Z}$ can be generated by the transformation $\mb {Z}=\mb {LW}$. Here is where the assumption of a Gaussian SRF is crucial. When $\mb {W} \sim N_ k(\Strong{0},\mb {I})$, then $\mb {Z}=\mb {LW}$ is also Gaussian. The mean of $\mb {Z}$ is

\[  E(\mb {Z}) = \mb {L}(E(\mb {W})) = \Strong{0}  \]

and the variance is

\[  \mbox{Var}(\mb {Z}) = \mbox{Var}(\mb {LW}) = E(\mb {LWW}’\mb {L}’) = \mb {L}E(\mb {WW}’)\mb {L}’ = \mb {LL}’ = \mb {C}  \]

Consider now an SRF $Z(\bm {s}), \bm {s}\in D \subset \mathcal{R}^2$, with spatial covariance function $C(\bm {h})$. Fix locations $\bm {s}_1,\bm {s}_2,\cdots ,\bm {s}_ k$, and let $\mb {Z}$ denote the random vector

\[  \mb {Z} = \left[ \begin{array}{c} Z(\bm {s}_1)\\ Z(\bm {s}_2)\\ \vdots \\ Z(\bm {s}_ k) \end{array} \right]  \]

with corresponding covariance matrix

\[  \mb {C}_ z = \left( \begin{array}{cccc} C(\Strong{0}) &  C(\bm {s}_1-\bm {s}_2) &  \cdots &  C(\bm {s}_1-\bm {s}_ k)\\ C(\bm {s}_2-\bm {s}_1) &  C(\Strong{0}) &  \cdots &  C(\bm {s}_2-\bm {s}_ k)\\ &  \ddots & & \\ C(\bm {s}_ k-\bm {s}_1) &  C(\bm {s}_ k-\bm {s}_2) &  \cdots &  C(\Strong{0}) \end{array} \right)  \]

Since this covariance matrix is symmetric and positive definite, it has a Cholesky root, and the $Z(\bm {s}_ i), i=1,\cdots ,k$, can be simulated as described previously. This is how the SIM2D procedure implements unconditional simulation in the zero-mean case. More generally,

\[  Z(\bm {s}) = \mu (\bm {s}) + \varepsilon (\bm {s})  \]

where $\mu (\bm {s})$ is a quadratic form in the coordinates $\bm {s}=(x,y)$ and the $\varepsilon (\bm {s})$ is an SRF that has the same covariance matrix $\mb {C}_ z$ as previously. In this case, the $\mu (\bm {s}_ i), i=1,\cdots ,k$, is computed once and added to the simulated vector $\varepsilon (\bm {s}_ i), i=1,\cdots ,k$, for each realization.

For a conditional simulation, this distribution of

\[  \mb {Z} = \left[ \begin{array}{c} Z(\bm {s}_1)\\ Z(\bm {s}_2)\\ \vdots \\ Z(\bm {s}_ k) \end{array} \right]  \]

must be conditioned on the observed data. The relevant general result concerning conditional distributions of multivariate normal random variables is the following. Let $\mb {X} \sim N_ m(\bmu ,\bSigma )$, where

\[  \mb {X} = \left[ \begin{array}{c} \mb {X}_1\\ \mb {X}_2\\ \end{array} \right]  \]
\[  \bmu = \left[ \begin{array}{c} \bmu _1\\ \bmu _2\\ \end{array} \right]  \]
\[  \bSigma = \left( \begin{array}{cc} \bSigma _{11} &  \bSigma _{12}\\ \bSigma _{21} &  \bSigma _{22}\\ \end{array} \right)  \]

The subvector $\mb {X}_1$ is $k\times 1$, $\mb {X}_2$ is $n\times 1$, $\bSigma _{11}$ is $k\times k$, $\bSigma _{22}$ is $n\times n$, and $\bSigma _{12}=\bSigma _{21}’$ is $k\times n$, with $k+n=m$. The full vector $\mb {X}$ is partitioned into two subvectors, $\mb {X}_1$ and $\mb {X}_2$, and $\bSigma $ is similarly partitioned into covariances and cross covariances.

With this notation, the distribution of $\mb {X}_1$ conditioned on $\mb {X}_2=\bm {x}_2$ is $N_ k(\tilde{\bmu },\tilde{\bSigma })$, with

\[  \tilde{\bmu }=\bmu _1+\bSigma _{12}\bSigma _{22}^{-1}(\bm {x}_2-\bmu _2)  \]


\[  \tilde{\bSigma } = \bSigma _{11}-\bSigma _{12}\bSigma _{22}^{-1}\bSigma _{21}  \]

See Searle (1971, pp. 46–47) for details. The correspondence with the conditional spatial simulation problem is as follows. Let the coordinates of the observed data points be denoted $\tilde{\bm {s}}_1,\tilde{\bm {s}}_2,\cdots , \tilde{\bm {s}}_ n$, with values $\tilde{z}_1,\tilde{z}_2,\cdots ,\tilde{z}_ n$. Let $\tilde{\mb {Z}}$ denote the random vector

\[  \tilde{\mb {Z}} = \left[ \begin{array}{c} Z(\tilde{\bm {s}}_1)\\ Z(\tilde{\bm {s}}_2)\\ \vdots \\ Z(\tilde{\bm {s}}_ n) \end{array} \right]  \]

The random vector $\tilde{\mb {Z}}$ corresponds to $\mb {X}_2$, while $\mb {Z}$ corresponds to $\mb {X}_1$. Then $\left(\mb {Z}\mid \tilde{\mb {Z}}=\tilde{\mb {z}}\right) \sim N_ k(\tilde{\bmu },\tilde{\mb {C}})$ as in the previous distribution. The matrix

\[  \tilde{\mb {C}} = \mb {C}_{11}-\mb {C}_{12}\mb {C}_{22}^{-1}\mb {C}_{21}  \]

is again positive definite, so a Cholesky factorization can be performed.

The dimension n for $\tilde{\mb {Z}}$ is simply the number of nonmissing observations for the VAR= variable; the values $\tilde{z}_1,\tilde{z}_2,\cdots ,\tilde{z}_ n$ are the values of this variable. The coordinates $\tilde{\bm {s}}_1,\tilde{\bm {s}}_2,\cdots ,\tilde{\bm {s}}_ n$ are also found in the DATA= data set, with the variables that correspond to the x and y coordinates identified in the COORDINATES statement. Note: All VAR= variables use the same set of conditioning coordinates; this fixes the matrix $\mb {C}_{22}$ for all simulations.

The dimension k for $\mb {Z}$ is the number of grid points specified in the GRID statement. Since there is a single GRID statement, this fixes the matrix $\mb {C}_{11}$ for all simulations. Similarly, $\mb {C}_{12}$ is fixed.

The Cholesky factorization $\tilde{\mb {C}}=\mb {LL}’$ is computed once, as is the mean correction

\[  \tilde{\bmu }=\bmu _1+\mb {C}_{12}\mb {C}_{22}^{-1}(\bm {x}_2-\bmu _2)  \]

The means $\bmu _1$ and $\bmu _2$ are computed using the grid coordinates $\bm {s}_1,\bm {s}_2,\cdots ,\bm {s}_ k$, the data coordinates $\tilde{\bm {s}}_1,\tilde{\bm {s}}_2,\cdots , \tilde{\bm {s}}_ n$, and the quadratic form specification from the MEAN statement. The simulation is now performed exactly as in the unconditional case. A $k \times 1$ vector of independent standard $N(0,1)$ random variables is generated and multiplied by $\mb {L}$, and $\tilde{\bmu }$ is added to the transformed vector. This is repeated N times, where N is the value specified for the NR= option.