Edgeworth Expansion Examples

Edgeworth expansions can be very tedious and complicated. The aim of this document is to illustrate them with concrete examples in the univariate setting and in the bivariate setting.

Univariate Edgeworth expansion

Suppose \(X_1,\dots,X_n\) are independent realizations of the random variable \(X\), where \(\mathbb{E} X = 0\), \(\mathbb{E}X^2 = 1\), \(\mathbb{E}X^3 =\gamma\), and \(\mathbb{E}X^4 = \tau < \infty\). Then \(\sqrt{n}\bar{X}_n\) converges in distribution to the \(\text{Normal}(0,1)\) distribution. The Edgeworth approximations \(f_n^{(1)}\) and \(f_n^{(2)}\) of to within \(o(n^{-1/2})\) and \(o(n^{-1})\), respectively, of the pdf and \(F_n^{(1)}\) and \(F_n^{(2)}\) to the cdf of \(\sqrt{n}\bar{X}_n\) are given by \[\begin{align*} f_n^{(1)}(x) &= \phi(x)\left[1 + \frac{\gamma}{6\sqrt{n}}H_3(x)\right]\\ f_n^{(2)}(x) &= \phi(x)\left[1 + \frac{\gamma}{6\sqrt{n}}H_3(x) + \frac{\tau - 3}{24 n}H_4(x) + \frac{\gamma^2}{72n}H_6(x)\right] \end{align*}\] and \[\begin{align*} F_n^{(1)}(x) &= \Phi(x) - \phi(x)\left[\frac{\gamma}{6\sqrt{n}}H_2(x)\right]\\ F_n^{(2)}(x) &= \Phi(x) - \phi(x)\left[\frac{\gamma}{6\sqrt{n}}H_2(x) + \frac{\tau - 3}{24 n}H_3(x) + \frac{\gamma^2}{72n}H_5(x)\right], \end{align*}\] where \(H_j(x)\), \(j = 0,\dots,6\) are the Hermite polynomials \[\begin{align*} H_0(x) &= 1\\ H_1(x) &= x\\ H_2(x) &= x^2 -1\\ H_3(x) &= x^3 - 3x\\ H_4(x) &= x^4 - 6x^2+3\\ H_5(x) &= x^5- 10x^3 + 15x\\ H_6(x) &= x^6 - 15 x^4 + 45 x^2 - 15 \end{align*}\] which are defined by the relation \[ (-1)^j\frac{d^j}{dx^j}\phi(x)=\phi(x)H_j(x) \] for all \(j\).

Example

Let \(Y_1,\dots,Y_n \overset{\text{ind}}{\sim}\text{Exponential}(\lambda)\). Then \(\sqrt{n}\bar{Y}_n \sim \text{Gamma}(n,\lambda/\sqrt{n})\). If we let \(X_i = (Y_i - \lambda)/\lambda\), then \(\sqrt{n}\bar{X}_n\) has density given by \[ f_n(x) = \frac{1}{\Gamma(n)(1/\sqrt{n})^n}(x+\sqrt{n})^{n-1}\exp\left(-\frac{x +\sqrt{n}}{1/\sqrt{n}}\right) \quad \text{ for } x > -\sqrt{n}. \] Moreover, \(\mathbb{E}X_1 = 0\), \(\mathbb{E}X_1^2 = 1\), and it can be shown that \(\gamma = \mathbb{E}X_1^3 = 2\) and \(\tau = \mathbb{E}X_1^4 = 9\). We can plug these quantities into the formulas above to get the Edgeworth expansions of orders \(2\) and \(3\). The following R code does this, and makes plots of the cdfs, the pdfs, and the approximation errors for both of them.

# define needed Hermite polynomials
H2 <- function(x){x^2 - 1}
H3 <- function(x){x^3 - 3*x}
H4 <- function(x){x^4 - 6*x^2 + 3}
H5 <- function(x){x^5 - 10*x^3 + 15*x}
H6 <- function(x){x^6 - 15*x^4 + 45*x^2 - 15}

# input values of gamma and tau
gamma <- 2
tau <- 9

# choose a sampe size
n <- 10

# get exact pdf and cdf
x <- seq(-4, 4, length = 500)
Exact_pdf <- dgamma(x + sqrt(n),shape = n, scale = 1/sqrt(n))
Exact_cdf <- pgamma(x + sqrt(n),shape = n, scale = 1/sqrt(n))

# get pdf and cdf of 2nd order Edgeworth approximation
EE1_pdf <- dnorm(x) * ( 1 + gamma / (6 * sqrt(n)) * H3(x))
EE1_cdf <- pnorm(x) - dnorm(x) * gamma / (6 * sqrt(n)) * H2(x)

# get pdf and cdf of 3rd order Edgeworth approximation
EE2_pdf <- EE1_pdf + dnorm(x) * ( (tau - 3) / (24*n) * H4(x) + gamma^2 / (72*n) * H6(x) )
EE2_cdf <- EE1_cdf - dnorm(x) * ( (tau - 3) / (24*n) * H3(x) + gamma^2 / (72*n) * H5(x) )

# get limiting cdf and pdf
Asymp_cdf <- pnorm(x)
Asymp_pdf <- dnorm(x)

# make plots
par(mfrow  = c(1,2))

# cdf approximation
plot(EE1_cdf ~ x, 
     type = "l",
     lty = 2,
     ylab = "cdf", 
     col = "red")

lines(EE2_cdf ~ x, lty = 2, col = "blue")
lines(Asymp_cdf ~ x, lty = 3)
lines(Exact_cdf ~ x)

# cdf approximation error
plot(Asymp_cdf  - Exact_cdf ~ x, 
     type = "l",
     lty = 3,
     ylab = "cdf approximation error")

abline(h = 0,lwd= .5)

lines(EE2_cdf  - Exact_cdf ~ x, lty = 2, col = "blue")
lines(EE1_cdf  - Exact_cdf ~ x, lty = 2, col = "red")

par(mfrow  = c(1,2))

# pdf approximation
plot(EE1_pdf ~ x, 
     type = "l",
     lty = 2,
     ylab = "cdf",
     ylim = range(EE1_pdf,EE2_pdf),
     col = "red")

lines(EE2_pdf ~ x, lty = 2, col = "blue" )
lines(Asymp_pdf ~ x, lty = 3)
lines(Exact_pdf ~ x)

# pdf approximation error
plot(Asymp_pdf  - Exact_pdf ~ x, 
     type = "l",
     lty = 3,
     ylab = "pdf approximation error")

abline(h = 0,lwd= .5)

lines(EE2_pdf - Exact_pdf ~ x, lty = 2, col = "blue")
lines(EE1_pdf - Exact_pdf ~ x, lty = 2, col = "red")

Edgeworth expansion of studentized mean

We now consider the distribution of \(\sqrt{n}\bar{X}_n/\hat{\sigma}_n\), where \(\hat{\sigma}^2_n =n^{-1}\sum_{i=1}^n(X_i - \bar{X}_n)^2\) when \(\mathbb{E} X = 0\), \(\mathbb{E}X^2 = 1\), \(\mathbb{E}X^3 =\gamma\), and \(\mathbb{E}X^4 = \tau < \infty\). In Peter Hall’s book we can find approximations to the pdf and cdf of \(\sqrt{n}\bar{X}_n/\hat{\sigma}_n\). We have \[\begin{align*} \hat f_n^{(1)}(x)& = \phi(x) \left[1 - \frac{\gamma}{6\sqrt{n}}(2x^3 - 3x)\right]\\ \hat f_n^{(2)}(x)& = \phi(x) \left[1 - \frac{\gamma}{6\sqrt{n}}(2x^3 - 3x) - \frac{\tau - 3}{12 n}(x^4 - 6x^2 + 3) + \frac{\gamma^2}{18 n}(x^6 - 3x^4 - 9 x^2 + 3) + \frac{1}{4 n}(x^4 - 3)\right]\\ \end{align*}\] and \[\begin{align*} \hat F_n^{(1)}(x) & = \Phi(x) + \phi(x) \left[\frac{\gamma}{6\sqrt{n}}(2x^2 + 1)\right]\\ \hat F_n^{(2)}(x) & = \Phi(x) + \phi(x) \left[\frac{\gamma}{6\sqrt{n}}(2x^2 + 1) + \frac{\tau - 3}{12 n}(x^3 - 3x) - \frac{\gamma^2}{18 n}(x^5 + 2x^3 - 3x) - \frac{1}{4 n}(x^3 + 3x)\right]. \end{align*}\]

To illustrate these approximations, we run a simulation in which we treat the women’s finishing times for the 2009 Boston Marathon as a population and repeatedly draw random samples of size \(n\). The data can be accessed with a url as in the code below:

load(url("https://people.stat.sc.edu/gregorkb/data/bos09whrs.Rdata"))

# compute moments
N <- length(bos09whrs)
mu <- mean(bos09whrs)
sigma <- sqrt(mean((bos09whrs - mu)^2))
gamma <- mean( (bos09whrs - mu)^3) /sigma^3
tau <- mean( (bos09whrs - mu)^4) /sigma^4

# make histogram
hist(bos09whrs,main = "",xlab = "Women's finishing times for 2009 Boston Marathon")

# choose sample size and number of simulations
n <- 15
S <- 100000

There are \(9303\) finishing times, the empirical distribution of which has \(\gamma = 1.09\) and \(\tau = 4.41\). Under \(n = 15\), we obtain an approximation to the true pdf and cdf of \(\sqrt{n}\bar{X}_n/\hat{\sigma}_n\) with a kernel density estimator based on \(10^{5}\) simulated samples. Then we compare the kernel density estimate of the pdf and cdf to the standard Normal pdf and cdf as well as to the two Edgeworth approximations.

X <- matrix(bos09whrs[sample(N,S*n,replace=TRUE)],S,n)
R <- sqrt(n)*(apply(X,1,mean) - mu)/apply(X,1,sd)

kde <- density(R,from = -4,to = 4)
kde_pdf <- kde$y
kde_cdf <- c(0,cumsum(kde$y[-1]*diff(kde$x)))
x <- kde$x

# get pdf and cdf of 2nd order Edgeworth approximation
EE1_pdf <- dnorm(x) * ( 1 - gamma / (6 * sqrt(n)) * (2*x^3 - 3*x))
EE1_cdf <- pnorm(x) + dnorm(x) * gamma / (6 * sqrt(n)) * (2*x^2 + 1)

# get pdf and cdf of 3rd order Edgeworth approximation
EE2_pdf <- EE1_pdf - dnorm(x) * ( (tau - 3) / (12*n) * (x^4 - 6*x^2 + 3) - gamma^2 / (18*n) * (x^6 - 3*x^4 - 9*x^2 + 3) - 1/(4*n)*(x^4 - 3) )
EE2_cdf <- EE1_cdf + dnorm(x) * ( (tau - 3) / (12*n) * (x^3 - 3*x) - gamma^2 / (18*n) * (x^5 + 2*x^3 - 3*x) - 1/(4*n)*(x^3 + 3*x) )

# get asymptotic pdf and cdf
Asymp_pdf <- dnorm(x)
Asymp_cdf <- pnorm(x)

# make plots
par(mfrow  = c(1,2))

# cdf approximation
plot(EE1_cdf ~ x, 
     type = "l",
     lty = 2,
     ylab = "cdf", 
     col = "red")

lines(EE2_cdf ~ x, lty = 2, col = "blue")
lines(Asymp_cdf ~ x, lty = 3)
lines(kde_cdf ~ x)

# cdf approximation error
plot(Asymp_cdf  - kde_cdf ~ x, 
     type = "l",
     lty = 3,
     ylab = "cdf approximation error",
     ylim = range(EE2_cdf  - kde_cdf,EE1_cdf  - kde_cdf, Asymp_cdf  - kde_cdf))

abline(h = 0,lwd= .5)

lines(EE2_cdf  - kde_cdf ~ x, lty = 2, col = "blue")
lines(EE1_cdf  - kde_cdf ~ x, lty = 2, col = "red")

# make plots
par(mfrow  = c(1,2))

# pdf approximation
plot(EE1_pdf ~ x, 
     type = "l",
     lty = 2,
     ylab = "cdf",
     ylim = range(EE1_pdf,EE2_pdf),
     col = "red")

lines(EE2_pdf ~ x, lty = 2, col = "blue" )
lines(Asymp_pdf ~ x, lty = 3)
lines(kde_pdf ~ x)

# pdf approximation error
plot(Asymp_pdf  - kde_pdf ~ x, 
     type = "l",
     lty = 3,
     ylab = "pdf approximation error")

abline(h = 0,lwd= .5)

lines(EE2_pdf - kde_pdf ~ x, lty = 2, col = "blue")
lines(EE1_pdf - kde_pdf ~ x, lty = 2, col = "red")

Multivariate Edgeworth expansion

Suppose \(X_1,\dots,X_n\) are independent realizations of the random variable \(X \in \mathbb{R}^d\), where \(\mathbb{E}X = 0\), \(\mathbb{E}XX^T = \boldsymbol{\Sigma}\). We will present only the Edgeworth expansion of order 2, which will depend on moments of order \(3\). Precisely, it depends on all of the moments \(\mathbb{E}[X_1^{\alpha_1}\cdots X_d^{\alpha_d}]\), where \(\alpha_j\) is a non-negative integer for \(j=1,\dots,d\) and \(\alpha_1+\dots+\alpha_d = 3\).

To give an expression for the multivariate Edgeworth expansion, it will be easier to make use of multi-indices: Let \(X^{\alpha} = X_1^{\alpha_1}\cdots X_d^{\alpha_d}\), and let \(|\alpha| =\alpha_1+\dots+\alpha_d\). Moreover, let \[ \partial^{\alpha}= \frac{\partial^{|\alpha|}}{\partial x_1^{\alpha_1}\cdots \partial x_d^{\alpha_d}}. \] Then the Edgeworth approximation to within \(o(n^{-1/2})\) of the joint density of \(\sqrt{n}\bar{X}_n\) is given by \[ f_n^{(1)}(x) = \phi(x;\boldsymbol{\Sigma})\left[1 + \frac{1}{6\sqrt{n}}\sum_{|\alpha|=3}\mathbb{E} X^\alpha \cdot H_\alpha(x;\boldsymbol{\Sigma})\right], \] where \(H_\alpha(x;\boldsymbol{\Sigma})\) are the multivariate Hermite polynomials defined by the relation \[ (-1)^{|\alpha|}\partial^{\alpha}\phi(x;\Sigma) = H_{\alpha}(x;\Sigma)\phi(x;\Sigma), \] where \(\phi(x;\boldsymbol{\Sigma})\) is the pdf of the multivariate Normal distribution with mean vector \(\mathbf{0}\) and covariance matrix \(\boldsymbol{\Sigma}\).

The bivariate case

When \(d=2\), the Edgeworth expansion involves moments of orders \(\alpha = (3,0),(0,3),(2,1),(1,2)\), that is \(\mathbb{E}X_1^3\), \(\mathbb{E}X_2^3\), \(\mathbb{E}X_1^2X_2\), and \(\mathbb{E}X_1X_2^2\), and the required Hermite polynomials are given by \[\begin{align*} H_{(3,0)}(x_1,x_2;\boldsymbol{\Sigma}) &= \Omega_{11}x_1^3 + 3\Omega_{11}^2\Omega_{12}x_2x_1^2 + 3\Omega_{11}\Omega_{12}^2x_1x_2^2 + \Omega_{12}^3 x_2^3 - 3\Omega_{11}^2 x_1 - 3\Omega_{12}\Omega_{11}x_2\\ H_{(0,3)}(x_1,x_2;\boldsymbol{\Sigma}) &= \Omega_{22}x_2^3 + 3\Omega_{22}^2\Omega_{21}x_1x_2^2 + 3\Omega_{22}\Omega_{21}^2x_2x_1^2 + \Omega_{21}^3 x_1^3 - 3\Omega_{22}^2 x_2 - 3\Omega_{21}\Omega_{22}x_1\\ H_{(2,1)}(x_1,x_2;\boldsymbol{\Sigma}) &= -(\Omega_{22}\Omega_{11} + 2\Omega_{12}^2)x_2 - 3\Omega_{11}\Omega_{12}x_1 + (\Omega_{22}\Omega_{11}^{2} + 2 \Omega_{11}\Omega_{12}^2)x_1^2x_2\\ & \quad \quad + (2\Omega_{11}\Omega_{12}\Omega_{22} + \Omega_{12}^3)x_1x_2^2 + \Omega_{22}\Omega_{12}^2x_2^3 + \Omega_{11}^2\Omega_{12}x_1^3\\ H_{(1,2)}(x_1,x_2;\boldsymbol{\Sigma}) &= -(\Omega_{11}\Omega_{22} + 2\Omega_{21}^2)x_1 - 3\Omega_{22}\Omega_{21}x_2 + (\Omega_{11}\Omega_{22}^{2} + 2 \Omega_{22}\Omega_{21}^2)x_2^2x_1 \\ & \quad \quad + (2\Omega_{22}\Omega_{21}\Omega_{11} + \Omega_{21}^3)x_2x_1^2 + \Omega_{11}\Omega_{21}^2x_1^3 + \Omega_{22}^2\Omega_{21}x_2^3 \end{align*}\] where \(\boldsymbol{\Omega}= \boldsymbol{\Sigma}^{-1}\).

Example

Let \(U_1,U_2,U_3 \overset{\text{ind}}{\sim}\text{Exponential}(\lambda)\) and let \(X = (X_1,X_2)^T\), where \(X_1 = U_1 + U_2 - 2\lambda\) and \(X_2 = U_1 + U_3 - 2\lambda\). Then \[ \text{Cov}(X) = \lambda^2\left[\begin{array}{cc} 2&1 \\ 1&2 \end{array}\right] =:\boldsymbol{\Sigma}, \quad \mathbb{E}X_1^3 = \mathbb{E}X_2^3 = 4\lambda^3,\quad \mathbb{E} X_1^2X_2 = \mathbb{E} X_1X_2^2 = 2\lambda^3. \] If \(X_1,\dots,X_n\) are independent realizations of \(X\), the Edgeworth approximation to the density of \(\sqrt{n}\bar{X}_n\) to within \(o(n^{-1/2})\) is given by \[\begin{align*} f_n^{(1)}(x_1,x_2) = \phi(x_1,x_2;\boldsymbol{\Sigma})&\left[1 + \frac{1}{6\sqrt{n}}\left( 4\lambda^3 [H_{(3,0)}(x_1,x_2;\boldsymbol{\Sigma}) + H_{(0,3)}(x_1,x_2;\boldsymbol{\Sigma})] \right.\right. \\ & \hspace{1in}\left.\left.+ 2\lambda^3 [H_{(2,1)}(x_1,x_2;\boldsymbol{\Sigma}) + H_{(1,2)}(x_1,x_2;\boldsymbol{\Sigma})]\right)\right]. \end{align*}\] The following R code computes a kernel density estimate of the true bivariate density \(f_n(x_1,x_2)\) of \(\sqrt{n}\bar{X}_n\) based on \(100000\) simulated data points and compares it to its Edgeworth approximation \(f_n^{(1)}(x_1,x_2)\) and to its limit density \(\phi(x_1,x_2;\boldsymbol{\Sigma})\). We compare the heights of the bivariate densities over a set of points \((x_1(t),x_2(t)) = (4t\cos(t4\pi),4t\sin(t4\pi))\) for \(t\in[0,1]\), which form a spiral going out from the origin. The heights of the bivariate densities over \(t\in[0,1]\) are plotted as well as the differences between the kernel density estimate of the density \(\sqrt{n}\bar{X}_n\) and the limiting density and the Edgeworth approximation.

library(latex2exp)

# necessary hermite polynomials
H_30 <- function(x,Omega){
        
        v1 <- Omega[1,2]^3 * x[2]^3
        v2 <- 3 * Omega[1,1]*Omega[1,2]^2 * x[1]*x[2]^2
        v3 <- 3 * Omega[1,2]*Omega[1,1]^2 * x[1]^2*x[2]
        v4 <- Omega[1,1]^3 * x[1]^3
        v5 <- - 3 * Omega[1,1]^2 * x[1]
        v6 <- - 3 * Omega[1,2]*Omega[1,1] * x[2]
        
        val <- v1 + v2 + v3 + v4 + v5 + v6
        
        return(val)
}

H_03 <- function(x,Omega){
        
        v1 <- Omega[2,1]^3 * x[1]^3
        v2 <- 3 * Omega[2,2]*Omega[2,1]^2 * x[2]*x[1]^2
        v3 <- 3 * Omega[2,1]*Omega[2,2]^2 * x[2]^2*x[1]
        v4 <- Omega[2,2]^3 * x[2]^3
        v5 <- - 3 * Omega[2,2]^2 * x[2]
        v6 <- - 3 * Omega[2,1]*Omega[2,2] * x[1]
        
        val <- v1 + v2 + v3 + v4 + v5 + v6
        
        return(val)
}

H_21 <- function(x,Omega){
     
        v1 <- -(Omega[2,2]*Omega[1,1] + 2*Omega[1,2]^2)  * x[2]
        v2 <- - 3*Omega[1,1]*Omega[1,2]  * x[1]
        v3 <- (Omega[2,2]*Omega[1,1]^2 + 2*Omega[1,1]*Omega[1,2]^2) * x[1]^2 * x[2]
        v4 <- (2*Omega[1,1]*Omega[1,2]*Omega[2,2] + Omega[1,2]^3 ) * x[1] * x[2]^2
        v5 <- Omega[2,2]*Omega[1,2]^2 * x[2]^3
        v6 <- Omega[1,1]^2*Omega[1,2] * x[1]^3
        
        val <- v1 + v2 + v3 + v4 + v5 + v6
        
        return(val)
}

H_12 <- function(x,Omega){
        
        v1 <- -(Omega[1,1]*Omega[2,2] + 2*Omega[2,1]^2)  * x[1]
        v2 <- - 3*Omega[2,2]*Omega[2,1]  * x[2]
        v3 <- (Omega[1,1]*Omega[2,2]^2 + 2*Omega[2,2]*Omega[2,1]^2) * x[2]^2 * x[1]
        v4 <- (2*Omega[2,2]*Omega[2,1]*Omega[1,1] + Omega[2,1]^3 ) * x[2] * x[1]^2
        v5 <- Omega[1,1]*Omega[2,1]^2 * x[1]^3
        v6 <- Omega[2,2]^2*Omega[2,1] * x[2]^3
         
        val <- v1 + v2 + v3 + v4 + v5 + v6
        
        return(val)
}

# bivariate limiting pdf
phi <- function(x,Sigma){
        
        val <- 1/(2*pi) * 1 / sqrt(det(Sigma)) * exp( - (1/2) * t(x) %*% solve(Sigma) %*% x)
        
        return(val)
}

# bivariate Edgeworth expansion
EE1 <- function(x,Sigma,n,E_30,E_03,E_21,E_12){
        
        Omega <- solve(Sigma)
        corr1 <- E_30 * H_30(x,Omega) + E_03 * H_03(x,Omega) + E_21 * H_21(x,Omega) + E_12 * H_12(x,Omega)
        val <- phi(x, Sigma)*( 1 + 1/(6*sqrt(n)) * corr1 )
        
        return(val)
}

# make a kernel density estimator:
biv_kde <- function(x,data,h){
        
        val <- mean(dnorm(data[,1] - x[1],0,h) * dnorm(data[,2]- x[2],0,h))
        return(val)
        
}

# now for a specific example:
lambda <- 1
n <- 10
E_30 <- 4*lambda^3
E_03 <- 4*lambda^3
E_12 <- 2*lambda^3
E_21 <- 2*lambda^3
Sigma <- lambda^2 * matrix(c(2,1,1,2),nrow = 2, byrow = TRUE)

# make spiral path
t.seq <- seq(0,1,length = 200)
xy.seq <- cbind(4*t.seq*cos(t.seq*4*pi),4*t.seq*sin(t.seq*4*pi))

# get Edgeworth approximation at points along spiral path
EE1.seq <- apply(xy.seq,1,FUN=EE1, Sigma = Sigma, n = n, E_30 = E_30, E_03 = E_03, E_21 = E_21, E_12 = E_12)

# get asymptotic density at points along spiral path
phi.seq <- apply(xy.seq,1,phi,Sigma=Sigma)

# simulate data
S <- 100000
U.bar <- matrix(rgamma(S*3,shape = n, scale = lambda/n),S,3)
X.bar <- sqrt(n) * cbind( U.bar[,1] + U.bar[,2] - 2*lambda, U.bar[,1] + U.bar[,3] - 2*lambda)

# get bivariate kde at points along spiral path
kde.seq <- apply(xy.seq,1,FUN = biv_kde,data = X.bar, h =.25)

# make contour plot of limiting density with spiral path overaid
x.seq <- seq(-4,4,length = 100)
y.seq <- seq(-4,4,length = 100)
z.mat <- matrix(0,length(x.seq),length(y.seq))
for( j in 1:length(x.seq))
        for(k in 1:length(y.seq)){
                
                z.mat[j,k] <- phi(c(x.seq[j],y.seq[k]),Sigma = Sigma)                
                
        }

contour(x = x.seq, 
        y = y.seq, 
        z = z.mat,
        nlevels= 30, drawlabels=FALSE,
        lty = 1,
        lwd=.5,
        main = TeX("Contours of $\\phi(x_1,x_2;\\Sigma)$ with path $(x_1(t),x_2(t)), 0 \\leq t \\leq 1$ overlaid"))
lines(xy.seq[,2] ~ xy.seq[,1],
      type = "l",
      col = "green",
      lwd = 2)

# plot joint densities over spiral path
plot(phi.seq ~ t.seq,
     ylim = range(EE1.seq,kde.seq,phi.seq),
     type = "l",
     xlab = TeX("$t$"),
     ylab = TeX("Height of joint density at $(x_1(t),x_2(t))$"),
     lty = 3)
lines(EE1.seq ~ t.seq, 
      col = "red",
      lty = 2)
lines(kde.seq ~ t.seq)

# plot approximation errors over spiral path
plot(phi.seq - kde.seq ~ t.seq,
     ylim = range(EE1.seq-kde.seq,phi.seq - kde.seq),
     type = "l",
     lty = 3,
     xlab = TeX("$t$"),
     ylab = TeX("Approximation error at $(x_1(t),x_2(t))$"))
lines(EE1.seq - kde.seq ~ t.seq,
      col="red",
      lty = 2)
abline(h = 0, lwd = .5)