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)$

and

$\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.