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,\ldots ,\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,\ldots ,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,\ldots ,k$, is computed once and added to the simulated vector $\varepsilon (\bm {s}_ i), i=1,\ldots ,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,\ldots , \tilde{\bm {s}}_ n$, with values $\tilde{z}_1,\tilde{z}_2,\ldots ,\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,\ldots ,\tilde{z}_ n$ are the values of this variable. The coordinates $\tilde{\bm {s}}_1,\tilde{\bm {s}}_2,\ldots ,\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,\ldots ,\bm {s}_ k$, the data coordinates $\tilde{\bm {s}}_1,\tilde{\bm {s}}_2,\ldots , \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.