lg: An R package for Local Gaussian Approximations

Abstract:

The package lg for the R programming language provides implementations of recent methodological advances on applications of the local Gaussian correlation. This includes the estimation of the local Gaussian correlation itself, multivariate density estimation, conditional density estimation, various tests for independence and conditional independence, as well as a graphical module for creating dependence maps. This paper describes the lg package, its principles, and its practical use.

Cite PDF Tweet

Published

Sept. 20, 2021

Received

Apr 4, 2020

Citation

Otneim, 2021

Volume

Pages

13/2

38 - 56


1 Introduction

propose the local Gaussian correlation (LGC) as a new measure of statistical dependence between two stochastic variables X1 and X2, which has the following important property yet unrivaled in the literature: It can separate between positive and negative nonlinear dependence while still reducing to the ordinary Pearson correlation coefficient if X1 and X2 are jointly normally distributed. The R-package localgauss provides two important functions in this context; one that calculates the sample LGC based on observed values of (X1,X2), and one that uses the estimated LGC to perform a local test of independence between X1 and X2 as described in detail by .

We have lately seen a number of new applications of the LGC that the localgauss package does not support, however. use the LGC to test for financial contagion across markets during crises. present a procedure for estimating multivariate density functions via the LGC, which modify in order to compute estimates of conditional density functions. present a test for serial independence within a time series, which extend in order to include a test for cross-correlation between two time series. Finally, develop the local Gaussian partial correlation (LGPC) as a measure of conditional dependence and a corresponding test for conditional independence.

This paper describes the lg package , which provides a unified framework to implement all these methods, as well as a tool for visualizing the LGC and LGPC as dependence maps. use the LGC in spectral analysis of time series, but those methods have their own computational ecosystem in the localgaussSpec package .

In Section 2, we provide a brief introduction to the LGC as well as the methods and applications referred to above. In Section 3, we describe the core function in the lg package and move on to demonstrate the implementation of various applications in Section 4. We conclude this paper in Section 5 by demonstrating the graphical capabilities of the lg package.

2 Statistical background

Consider a random vector X having the unknown probability density function fX(x). It is a standard task to estimate fX based on a random sample X1,,Xn, and the statistical literature provides an abundance of methods to accomplish this. One may, for example, make the assumption that the unknown density function has a particular parametric form, fXFθ, where Fθ={f(x;θ),θΘ} is a family of probability density functions indexed by some parameter θ, and where Θ is the parameter space. Under this assumption, we will typically produce an estimate of the parameter θ, written θ^, using the maximum likelihood method. The estimated probability density function is then given as f^X(x)=f(x;θ^).

A different approach is to estimate fX() without any prior parametric assumptions. The classical method for nonparametric density estimation is the kernel estimator f^X(x)=1nbi=1nK(Xixb), where K is a symmetric density function (the kernel) and b is a tuning parameter (the bandwidth) that controls the smoothness of the estimate f^X(). See for an introduction to this topic. There is also a massive statistical literature on density estimation containing extensions and improvements to the classical methods to be used in various practical situations.

provide one such idea. They consider a parametric family Fθ, but instead of searching for a single parameter value θ0 for which fX(x)=f(x;θ0) (or approximately so), they rather assert that different members of Fθ may approximate fX locally in different parts of its domain. In other words, they seek to estimate a parameter function θ0(x) for which fX(x)=f(x;θ0(x)) (or approximately so), and do this by maximizing a local likelihood function in each point x: (1)θ^(x)=argmaxθΘLn(θ,x)=argmaxθΘ1nbi=1nK(Xixb)logf(Xi;θ)1bK(yxb)f(y;θ)dy, where, again, K is a symmetric density function and b is a bandwidth parameter that controls the smoothness of the estimate. The second term in the local likelihood function is a penalty that ensures that the estimated density f^X(x)=f(x;θ^(x)) converges correctly to the true density function fX(x) as the sample size n increases to infinity and the bandwidth b decreases towards zero. See for a detailed discussion about this construction.

consider the bivariate case X=(X1,X2) and take Fθ to be the family of bivariate normal distributions consisting of densities on the form

(2)f(x;θ)=ψ(x1,x2;μ1,μ2,σ1,σ2,ρ)=12πσ1σ21ρ2×exp{12(1ρ2)((x1μ1)2σ122ρ(x1μ1)(x2μ2)σ1σ2+(x2μ2)2σ22)}, where θ=(μ1,μ2,σ1,σ2,ρ) is the vector of parameters. Using a sample {X1i,X2i}, i=1,,n, they estimate θ locally in the point x by maximizing the local likelihood function (1), producing θ^(x)=(μ^1(x),μ^2(x),σ^1(x),σ^2(x),ρ^(x)), and take special interest in the estimated correlation function ρ^(x) (i.e., the LGC) because it serves as an attractive local measure of statistical dependence between X1 and X2. They show that the LGC reveals many types of nonlinear statistical dependence that are not captured by the ordinary (global) Pearson correlation coefficient. Furthermore, the LGC distinguishes between positive and negative dependence and reduces to the Pearson ρ if X1 and X2 are jointly normal. We refer to for a detailed treatment of the theoretical foundations of the LGC as well as several examples and rather present two simple illustrations at this point in order to demonstrate the concept.

graphic without alt text graphic without alt text
  1. Estimated local Gaussian correlations between two jointly normal variables having correlation equal to 0.5, based on 1000 observations.
  1. Local Gaussian correlation between daily returns from the CAC40 and FTSE100 stock indices based on 1000 consecutive trading days.
Figure 1: Two dependence maps

In Figure 1, we have plotted the estimated LGC for two bivariate data sets on a grid; 1000 simulated observations from a binormal distribution having correlation equal to 0.5, and the daily return on the CAC40 and FTSE100 stock indices on 1000 consecutive trading days starting on May 5th 2014 . In the first panel, we see that the estimated local correlation coincides with the global correlation, except for the estimation error which is comparable to the uncertainty observed in other nonparametric estimation methods such as the kernel density estimator (see, for instance, for a formal asymptotic analysis of relevant convergence rates). In the second panel, we see clearly that the local correlation, and thus the dependence, is stronger in the lower left and upper right regions of the distribution than in the central parts. The phenomenon of local dependence is well known in the financial literature, and using the LGC it can be measured, interpreted, and visualized in a natural way. The interpretation of this particular figure is that extreme observations on the two stock indices are more strongly dependent than the less extreme observations.

One may obtain these particular estimates from the older localgauss package (as well as the lg package, of course), but the plotting routine that was used to produce these figures is included in the lg package and will be described in more detail in Section 5.

Taking the LGC as a measure of dependence opens up a number of possibilities to construct statistical tests. show that ρ(x)0 implies that X1 and X2 are independent. They show further that independence between X1 and X2 implies ρ(x)0 if the population values of the local mean and standard deviation functions satisfy the following conditions: μi(x1,x2)=μi(xi) and σi(x1,x2)=σi(xi) for i=1,2. Equivalence between independence and ρ(x)0 holds in general if the observations have been suitably transformed according to a procedure presented later in this section. It follows then that departures from ρ^(x)0 may be taken as evidence against the hypothesis that X1 and X2 are statistically independent. formalize this notion by testing whether ρ(x)0 for all xSR2 using the test statistic

(3)Tn,b=Sh(ρ^(x))dFn(x) for some non-negative function h, for example h(x)=x2 or h(x)=|x|. Critical values may be obtained by permutations of the data under the null hypothesis, and we demonstrate the implementation of this test using the lg package in Section 4.2.

Consider next the stationary time series {Xt}. The autocorrelation function (ACF) ρk=ρ(Xt,Xtk) is a well known concept for describing the serial dependence in the time series, but the ACF is, again, only capable to completely capture linear serial dependence. seek to remedy this by rather calculating the local correlation between Xt and Xtk. This leads to a test for serial independence in a natural way. In fact, this work is mainly a theoretical exercise in order to accommodate time series dependence. Testing for independence between Xt and Xtk using observations {Xt,Xtk}t=k+1T leads to the same test statistic (3) and bootstrap procedure as the test for independence between X1 and X2 that we described above.

extend this problem to test for serial cross-dependence between two time series {Xt,Yt} by measuring the LGC between Xt and Ytk. Departures from ρ^(x,y)0 are, again, taken as evidence against independence, and the test statistic ((3)) provides an aggregate measure of this discrepancy in the specified region S. In this case, however, we can not obtain replicates of the test statistic under the null hypothesis by simple permutations of the data. suggest two block bootstrap procedures instead to this end, using fixed and random block sizes, respectively. The tests for serial dependence and serial cross-dependence are both implemented in the lg package, as we demonstrate in Section 4.2.

We find another application of the local Gaussian approximation in work by , who measure and test for financial contagion. They define contagion as "a significant increase in cross-market linkages after a shock to one country" and employ the LGC to quantify this potential linkage. The authors estimate the LGC on a grid {x1,x2}k, k=1,,K along the diagonal D={(x1,x2):x1=x2} before and after some critical event in the financial markets, denoted as the crisis (C) and the non-crisis (NC) periods, respectively. They compare the two estimates using the following test statistic, Tn,bD=k=1K{ρ^C(xk,xk)ρ^NC(xk,xk)}w(xk,xk), where w(,) is a weight function that serve the same purpose as the integration area S in ((3)). In this case, show that a standard bootstrap will suffice in order to produce approximate replicates of Tn,bD under the null hypothesis of no financial contagion, and we demonstrate the implementation of this test using the lg package in Section 4.3.

Although the original work by provide a general framework for local likelihood density estimation using any p-variate parametric family as the local family, it is evident that the method may struggle in multivariate applications much in the same way as the kernel density estimator does. This is a consequence of the curse of dimensionality, the effect of which is sought remedied by an algorithm provided by . The idea is to fit the p-variate normal distribution ψ(μ,Σ) locally, where μ is the vector of p expectations, and Σ is the p×p covariance matrix (to which the correlation matrix R corresponds), but under the following structural simplifications:

(4)μi(x)=μi(x1,,xp)=defμi(xi)

(5)σi(x)=σi(x1,,xp)=defσi(xi)

(6)ρij(x)=ρij(x1,,xp)=defρij(xi,xj). estimate the local parameters above by first obtaining univariate marginal locally Gaussian fits (eqs. (4) and (5)), and then pairwise bivariate locally Gaussian fits (eq. (6)). In the second step, the estimates μ^i(xi), μ^j(xj), σ^i(xi), and σ^j(xj) are kept fixed in the estimation of the pairwise local correlation. They argue further that the following transformation technique will produce better density estimates in many situations. The motivation for introducing the simplifications defined in equations (4)-(6) can be compared to the practical advantages of estimating additive regression models, where E(Y)=f(x1,,xp)=deff1(x1)++fp(xp).

Denote by Fi(xi), i=1,,p the marginal distribution functions of the stochastic vector X, and by F^i(xi)=n1i=1n1(Xixi) their empirical counterparts. They then estimate the density fZ(z) of the vector Z={Φ1(Fi(Xi))}i=1,,p. In practice it is approximated by

(7)Z^={Φ1(F^i(Xi))}i=1,,p, and where Ψ() is the univariate standard normal cdf. In that case, they simplify the estimation problem even further and fix

(8)μi(zi)=def0 and σi(zi)=def1,i=1,,p, so that the only parameter functions left to estimate are the pairwise local Gaussian correlations R(z)={ρij(zi,zj)}i<j. We use the notation Z, zi, and zj to signify that the estimation is performed on the (approximate) standard normal scale or z-scale for short. We can then estimate the joint density fZ(z) of Z as (9)f^Z(z)=ψ(z;μ(z)=0,σ(z)=1,R=R^(z)), where μ(z)={μi(z)} and σ(z)={σi(z)} for i=1,,p, and then substitute fZ for f^Z in the following relation obtained by in order to estimate the unknown density fX: (10)f(x)=fZ(Φ1(F1(x1)),,Φ1(Fp(xp)))×i=1pfi(xi)ϕ(Φ1(Fi(xi))), where ϕ() is the standard normal pdf. This estimator is implemented the lg package as demonstrated in Section 4.1.

One particular feature enjoyed by the jointly normally distributed vector X is that for any partitioning X=(X(1),X(2)), the conditional distribution of X(1)|X(2)=x(2) is also normal. In fact, if XN(μ,Σ), and μ and Σ is partitioned according to (X(1),X(2)) as μ=(μ1μ2) and Σ=(Σ11Σ12Σ21Σ22), then X(1)|X(2)=x(2)N(μ,Σ), where

(11)μ=μ1+Σ12Σ221(x(2)μ2)

(12)Σ=Σ11Σ12Σ221Σ21, see e.g. . demonstrate that this property may be translated into a corresponding local argument without modification. That is, if the joint density fX() can be written on a locally Gaussian form fX(x)=ψ(x,μ(x),Σ(x)), then the conditional density of X(1)|X(2)=x(2) can also be written on the same locally Gaussian form with local parameters given by equations ((11)) and ((12)), except that all quantities are x-dependent. If we use the transformation technique described above together with simplification ((8)), the local versions of equations ((11)) and ((12)) simplify to

(13)μ(z)=R12(z)R22(z)1z(2),

(14)Σ(z)=R11(z)R12(z)R22(z)1R21(z), where we, again, switch to z-notation in order to make it clear that these quantities are estimated on the standard normal z-scale. An estimator f^X(1)|X(2)(|) of the conditional density fX(1)|X(2)(|) follows immediately from an expression corresponding to ((10)), and the lg package provides functions for implementing this estimator in R. We describe the implementation of this functionality in Section 4.1.

Finally, we refer to who take the local version of the conditional covariance matrix ((12)) (or ((14)) in the transformed case) as a measure for conditional dependence, and thus as an instrument to test for conditional independence. Consider the stochastic vector X=(X(1),X(2),X(3)), where X(1) and X(2) are scalars and X(3) may be a vector. X(1) is conditionally independent from X(2) given X(3), written X(1)X(2)|X(3), if the stochastic variables X(1)|X(3) and X(2)|X(3) are independent, or, equivalently, if the joint conditional density of X(1) and X(2) given X(3) can be written as the product

(15)fX(1),X(2)|X(3)(x(1),x(2)|x(3))=fX(1)|X(3)(x(1)|x(3))×fX(2)|X(3)(x(2)|x(3)). In this case, denote by α(z) the off-diagonal element in the 2×2 local correlation matrix R(z) (which derives directly from Σ(z) as given in ((14))). If X has a local Gaussian distribution, the conditional independence ((15)) is equivalent to α(z)0, and take departures from this relation as evidence against the hypothesis of conditional independence between X(1) and X(2) given X(3). The natural way to quantify this is the test functional

(16)Tn,bCI=h(α^(z))dFn(z). describe a bootstrap procedure for generating replicates of Tn,bCI under the null hypothesis. In Section 4.4, we demonstrate how the lg package may be used to extract estimates of the local partial correlation and perform tests for conditional independence according to this scheme.

3 The first step: Creating the lg-object

The local Gaussian correlation may be used to perform a number of statistical analyses, as is evident from the preceding section. The practitioner must first, however, make three quite specific modeling choices; namely (i) to choose an estimation method, i.e., the level of simplification in multivariate applications, (ii) to determine whether the data should be transformed towards marginal standard normality before estimating the LGC, and (iii) to choose a set of bandwidths or at least a method for calculating bandwidths. The architecture of the lg package requires the user to make these choices before endeavoring further into specific applications by imposing a strict, two-step procedure:

  1. Create an lg-object.

  2. Apply relevant analysis functions to the lg-object.

In the following, we assume that one has a data set x loaded into the R workspace, which must be an n×p matrix (one column per variable, one row per observation), possibly including NAs which will be excluded from the analysis, or a data frame having the same dimensions. The fundamental syntax for creating an lg-object is lg_object <- lg_main(x), and we will, in this section, explain how the modeling decisions (i)-(iii) can be encoded into the lg-object by using appropriate arguments in this function.

Selecting the estimation method

Given a data set x having n rows and p2 columns, the user must choose between four distinct estimation methods and specify this choice by using the argument est_method to the lg_main()-function. We look at the built-in bivariate data set faithful, which records the waiting time between eruptions and the duration of the eruption for the Old Faithful geyser in the Yellowstone National Park, USA (see the help file in R for more details: ?faithful), and load the lg package in order to demonstrate the implementation:

R> x <- faithful
R> library(lg)

1. A full locally Gaussian fit for bivariate data. If p=2, we may fit the bivariate normal ψ(x,θ) locally to f(x), and by a "full local fit", we mean that we jointly estimate the five parameters θ(x)=(μ1(x1,x2),μ2(x1,x2),σ1(x1,x2),σ2(x1,x2),ρ(x1,x2)) by optimizing the local likelihood function ((1)) in the grid point x=(x1,x2). To use this estimation method in the subsequent analysis, specify est_method = "5par" in the call to lg_main():

R> lg_object <- lg_main(x, est_method = "5par")

The resulting lg_object is a list of class lg, and we may confirm that the assignment has been carried out correctly by inspecting its est_method-element:

R> lg_object$est_method 

[1] "5par"

Note that the full locally Gaussian fit for raw data is not available if the number of variables p is greater than 2. The lg_main()-function will check for this and print out an error message if p>2 and est_method = "5par".

2. A simplified locally Gaussian fit for multivariate data. As described in the preceding section, we may construct a simplified estimation procedure for calculating the LGC in two steps, which in principle works for any number of dimensions (including p=2):

  1. Calculate μi(xi) and σi(xi), i=1,,p by fitting the univariate normal distribution locally to each marginal density fi(xi) of f(x).

  2. Keep μ^i(x) and σ^i(xi), i=1,,p from step 1 fixed when estimating ρij(xi,xj), 1i<jp by fitting the bivariate normal distribution to each pair of variables Xi and Xj.

To use this method, create the lg-object by running the following line:

R> lg_object2 <- lg_main(x, est_method = "5par_marginals_fixed")

3. A simplified locally Gaussian fit for marginally standard normal data. This estimation method is applicable for marginally standard normal data, or data that have been transformed to approximate marginal standard normality by, e.g., the transformation ((7)). In that case, we fix the marginal expectation functions and standard deviation functions to the constant values 0 and 1, respectively, and estimate only the pairwise local Gaussian correlations as in ((9)). To use this estimation method, create the lg-object by running

R> lg_object3 <- lg_main(x, est_method = "1par")

Note that the function call above will issue a warning if the option for transforming the data to marginal standard normality is not at the same time set to TRUE, see the next sub-section on data transformation for details.

4. A full locally Gaussian fit for trivariate data. If the number of variables p is equal to 3, and we choose to transform the data to marginal standard normality (see the next sub-section), the transformed density fZ() in ((9)) may be estimated by jointly estimating the three local correlations ρ12(z1,z2,z3), ρ13(z1,z2,z3), and ρ23(z1,z2,z3). This estimation method was introduced recently by in order to increase power of their conditional independence test, but it can be used in any application described in this paper that consider trivariate data. To use this estimation method, create the lg-object by running

R> lg_object4 <- lg_main(x, est_method = "trivariate")

This command will throw an error if the data set x does not have exactly three columns.

Data transformation

graphic without alt text graphic without alt text
  1. The original data.
The same data transformed to marginal standard normality.
Figure 2: The same data on different scales.

Next, the user must determine if the local Gaussian correlation should be estimated directly on the raw data or on the marginally normal pseudo observations (7). This is carried out by using the logical transform_to_marginal_normality-argument in lg_main, for example:

R> lg_object <- lg_main(x, transform_to_marginal_normality = TRUE)

The resulting lg_object now includes the element transform_to_marginal_normality set according to the input, and if this is TRUE, it also includes the transformed_data and a function trans_new() that may be used later to apply the same transformation to, e.g., grid points. If the transformation option is set to FALSE, the transformed_data element contains the input data x, and trans_new() is nothing more than the identity mapping for points in Rp. See Figure 2 for two plots that demonstrate the effect of the data transformation on the example data.

Bandwidth selection

Finally, the user must specify a set of bandwidths or a method for calculating them. Given that the different estimation methods described in Section 3.1 require different sets of bandwidths (i.e, joint, marginal, and/or pairwise), the easiest approach for the user is to leave the selection and formatting of the bandwidths to the lg_main()-function.

The bandwidth plays a slightly different role in local likelihood estimation than elsewhere in the nonparametric literature. It controls the level of localization and thus only indirectly the smoothness of the estimates. Indeed, suppose we concentrate on the univariate case for the moment and assume that the (single) bandwidth b is small. In that case, we see from the local likelihood function (1) that only the very few observations closest to a fixed grid point x0 will have significant weight when determining the local parameter estimates θ^0(x) at that point. Moving on to another nearby point, x1 may then lead to a fairly different estimate θ^(x1) because the set of observations having weight in this point is very different. This may, again, lead to rougher parameter estimates θ^(x) and in turn also to rougher density estimates f(x,θ^(x)).

If the bandwidth b grows large, on the other hand, all observations receive similar weights, and furthermore: the local likelihood function (1) becomes approximately proportional to the ordinary (global) likelihood function Ln(θ)=i=1nlogf(Xi;θ). In other words, the local parameter estimates θ^(x) are smoothed towards the constant maximum likelihood estimates θ^ML, and the estimated density f(x;θ^(x)) towards the maximum likelihood estimate f(x;θ^ML). This means that the bandwidth may be chosen to reflect the goodness-of-fit of f(x;θ) to the true density f(x).

In the multivariate applications referred to in this paper, the bandwidth b in (1) is a diagonal matrix, and 1/b is naturally taken to represent its inverse.

We have in practice seen two automatic bandwidth selectors employed in the applications referred to in Section 2: a cross-validation procedure that is fairly slow to compute but accurate with respect to density estimation, and a plug-in bandwidth that is much quicker to calculate but less accurate with respect to density estimation. We use the argument bw_method to the lg_main()-function in order to choose between the two.

1. Choosing bandwidths by cross-validation. The functional

CV(b)=1ni=1nlogf(Xi;θ^(i)(Xi)), where θ^(i)(x) is the parameter estimate obtained after deleting observation Xi from the data, is proportional to a quantity that estimates the Kullback-Leibler distance between f(,θ^()) and the true density f(); see . The cross-validated bandwidth bCV is hence given by bCV=argminbCV(b). If we, for example, wish to use the simplified estimation procedure on the transformed data, we need bandwidths for the marginal estimates of the local means and standard deviations, as well as a 2×2 diagonal bandwidth matrix for each pair of variables. This is accomplished by the following call to lg_main():

R> # Create the lg-object with bandwidths chosen by cross-validation
R> lg_object <- lg_main(x,
R+                      est_method = "5par_marginals_fixed",
R+                      transform_to_marginal_normality = TRUE,
R+                      bw_method = "cv")

The lg_object now contains the necessary bandwidths for this configuration, as can be seen by inspecting the contents of its bw-element:

R> # Print out the bandwidths
R> lg_object$bw

$marginal
[1] 0.9989327 0.9875333

$marginal_convergence
[1] 0 0

$joint
x1 x2         bw1      bw2 convergence
1  1  2 0.2946889 0.331971           0

This is itself a list, containing the crucial elements marginal for the p marginal bandwidths, and joint that contains the p(p1)/2 bandwidth matrices, one for each pair of variables (which in this bivariate example just one variable pair, (x1, x2)). The convergence flags stem from the built-in R functions optim() and optimize() that we use to obtain the minimizer of CV(), and 0 indicates successful convergence.

2. Using plug-in bandwidths. Obtaining cross-validated bandwidths is unfortunately fairly slow on a standard computer. For sample sizes in the 500-1000 range, the process may take several minutes, which is unfeasible when embarking on analyses that require, e.g., resampling. We have, therefore, implemented a quick plug-in bandwidth selector as well that may suffice in many practical situations, especially at the initial or exploratory stage.

show that the simplified version of the local Gaussian fit have the same convergence rates as the corresponding nonparametric kernel density estimator for which derives the plug-in formula b=1.08sd(x)n1/5. By specifying bw_method = "plugin", the lg_main()-function will select the bandwidths correspondingly, except that the exponent changes to 1/6 for joint bandwidths, and the proportionality constant is by default set to 1.75. The latter number is the result of regressing bCV on n1/6 in a large simulation experiment covering various data generating processes . We see the effect of switching to plug-in bandwidths in the code below:

R> # Make the lg-object with plugin bandwidths
R> lg_object <- lg_main(x,
R+                      est_method = "5par_marginals_fixed",
R+                      transform_to_marginal_normality = TRUE,
R+                      bw_method = "plugin")
R> # Print out the bandwidths
R> lg_object$bw

$marginal
[1] 0.5703274 0.5703274

$marginal_convergence
[1] NA NA

$joint
  x1 x2       bw1       bw2 convergence
1  1  2 0.6875061 0.6875061          NA

Summary of the initialization function

Table 1: Arguments to the initialization function lg_main()
Argument Explanation Default value
x The data, an n×p matrix or data frame
bw_method Method for calculating the bandwidths "plugin"
est_method Estimation method "1par"
transform_to_
marginal_normality Transform the data TRUE
bw The bandwidths to use if already calculated NULL
plugin_constant_
marginal Prop. const. in plugin formula for marg. bw. 1.75
plugin_exponent_
marginal Exponent in plugin formula for marg. bw. 1/5
plugin_constant_
joint Prop. const. in plugin formula for joint bw. 1.75
plugin_exponent_
joint Exponent in plugin formula for joint bw. 1/6
tol_marginal Abs. tolerance when optimizing CV(b), marg. 103
tol_joint Abs. tolerance when optimizing CV(b), joint 103

In the sub-section above, we present the three most important arguments to lg_main(). Each of them allows the user to configure one of the three crucial modeling choices. Let us complete this treatment by covering some possibilities to make further adjustments to those choices.

  1. The user may supply the bandwidths directly to lg_main() by passing them to the bw-argument. They have to be in the correct format, though, which is a list containing the vector $marginal if est_method = "5par_marginals_fixed", and always a data frame $joint specifying all variable pairs in the x1 and x2 columns and the corresponding bandwidths in the bw1 and bw2 columns. The function bw_simple() will assist in creating bandwidth objects.

  2. If bw_method = "plugin" the user may change the proportionality constant and exponent in the plugin formula for the joint and, if applicable, the marginal bandwidths. See Table 1 for the necessary argument names.

  3. If bw_method = "cv", the user may change the numerical tolerance in the optimization of CV(b). See Table 1 for the necessary argument names.

4 Statistical inference using the lg package

We proceed in this section to demonstrate how to implement each of the tasks that we discussed in Section 2. The general pattern is to pass the lg-object to one of the estimation or test functions provided in the lg package. We will look at some financial data in the examples: the monthly returns on the S&P500, FTSE100, DAX30, and TOPIX stock indices from January 1985 to March 2018 .

Density estimation

We start by introducing a basic function for estimating the LGC on a grid as described by , and thus also a probability density estimate. We create a grid, x0, having the same number of columns as the data in the code below. Note that we use the pipe operator %>% from the magrittr package as well as functions from the dplyr package for easy manipulation of data frames. We then pass the grid and the lg-object containing our modeling choices to the dlg()-function in order to do the estimation.

graphic without alt text graphic without alt text
  1. The estimated pairwise LGC for the stock data along the diagonal indicating that the pairs of stocks are more strongly dependent in the lower part of the distribution than the upper part.
  1. The corresponding joint density estimate of the four return series, plotted along the diagonal in 4.
Figure 3: Local correlations and density estimates calculated using the dlg()-function.
R> # Create an lg-object
R> lg_object <- lg_main(x = stock_data %>% select(-Date),
R+                      est_method = "1par",
R+                      bw_method = "plugin",
R+                      transform_to_marginal_normality = TRUE)
R> 
R> # Construct a grid diagonally through the data.
R> grid_size <- 100
R> x0 <- stock_data %>% 
R+   select(-Date) %>% 
R+   apply(2, function(y) seq(from = -7,
R+                            to = 7,
R+                            length.out = grid_size))
R> 
R> # Estimate the local Gaussian correlation on the grid
R> density_object <- dlg(lg_object, grid = x0)

The last line of code creates a list containing a number of elements. The two most important are $loc_cor, which is a matrix of local correlations having one row per grid point and one column per pair of variables (the columns correspond to the rows in density_object$pairs), and $f_est, which is a vector containing the estimate f^X(x) of the joint density fX(x) in the grid points specified in x0. The estimated local correlations for this example is plotted in Figure 3a, and the corresponding density estimate is plotted (along the diagonal x1=x2=x3=x4=x) in Figure 3b.

The list density_object contains the estimated standard deviations of the local correlations in $loc_cor_sd, as well as lower and upper confidence bands $loc_cor_lower and $loc_cor_upper at the 95% level. We refer to Table 2 for a complete overview of the arguments to dlg().

Note that the configuration transform_to_marginal_normality = TRUE and est_method = 5par in the bivariate case coincides with the situation considered by . In that case, dlg() serves as a wrapper for the function localgauss() in the localgauss-package .

Obtaining the estimate of a conditional density using the algorithm described in Section 2 is very similar. However, one must take particular care of the ordering of the variables in the data set. The estimation function, clg(), will always assume that the free variables come first and the conditioning variables last. Let us illustrate this in the following code chunk by estimating the joint conditional density of S&P500 and FTSE100, given that DAX30 = TOPIX = 0.

Table 2: Arguments to the dlg()-function
Argument Explanation Default value
lg_object The lg-object created by lg_main()
grid The evaluation points for the LGC NULL
level Level for confidence bands 0.95
normalization The estimated density does not integrate to one by
_points c onstruction. dlg() will generate the given number
of normal variables, having the same moments as
the data, approximate f^X(x)dx by a Monte Carlo
integral, and then normalize the density estimate
accordingly
NULL
bootstrap Calculate bootstrapped confidence intervals instead
of asymptotic expressions FALSE
B Number of bootstrap replicates 500
R> # We must make sure that the free variables come first
R> returns1 <- stock_data %>% select(SP500, FTSE100, DAX30, TOPIX)
R> 
R> # Create the lg-object
R> lg_object <- lg_main(returns1,
R+                      est_method = "1par",
R+                      bw_method = "plugin",
R+                      transform_to_marginal_normality = TRUE)
R> 
R> # Create a grid
R> x0 <- returns1 %>% 
R+   select(SP500, FTSE100) %>% 
R+   apply(2, function(y) seq(from = -7,
R+                            to = 7,
R+                            length.out = grid_size))
R> 
R> # Calculate the conditional density
R> cond_density <- clg(lg_object, grid = x0, cond = c(0, 0))

The key argument in the call to clg() above is cond = c(0, 0). This means that the last two variables are conditioning variables (and hence, that the first 42=2 variables are free). The value of the conditioning variables are fixed at DAX30 = 0 and TOPIX = 0, respectively. This also means that the number of columns in the grid x0 plus the number of elements in cond must equal the number of variables p in the data set, and the call to clg() will result in an error message if this requirement is not fulfilled. The clg()-function takes mostly the same arguments as dlg() listed in Table 2, and the conditional density estimate in our example is available in the vector cond_density$f_est.

Tests for independence

Three independence tests based on the LGC have appeared in the literature thus far:

  1. A test for independence between the stochastic variables X1 and X2 based on iid data, cf. .

  2. A test for serial independence between Xt and Xtk within a time series {Xt}, cf. .

  3. A test for serial cross-dependence between Xt and Ytk for two time series {Xt} and {Yt}, cf. .

As we noted in Section 2, their practical implementations are very similar, and the lg package provides the function ind_test() to perform the tests. Let us first consider the i.i.d. case, and generate 500 observations test_x from the well known parabola model X2=X12+ε, where both X1 and ε are independent and standard normal. In this case, X1 and X2 are uncorrelated but obviously strongly dependent. considers mainly the full bivariate fit to the raw data, which we easily encode into the lg-object as before. The implementation of the test using 100 bootstrap replicates is shown below.

R> # Make the lg-object
R> lg_object <- lg_main(test_x, 
R+                      est_method = "5par",
R+                      transform_to_marginal_normality = TRUE)
R> # Perform the independence test
R> test_result <- ind_test(lg_object, n_rep = 100)
R> # Print out the p-value of the test
R> test_result$p_value

[1] 0

This may take a few minutes to run on a desktop computer due to bootstrapping. The small p-value indicates that we reject the null hypothesis of independence between X1 and X2 in the parabola model defined above. We can further specify the function h and the integration area S in the test statistic (3); see Table 3 for details.

Table 3: Arguments to the ind_test()-function
Argument Explanation Default value
lg_object The lg-object created by lg_main()
h The function h() in ((3)) function(x) x``^``2
S The integration area S in ((3)). Must function(x)
be a logical function on potential as.logical(rep(1,
grid points in R2 nrow(x)))
bootstrap The bootstrap method, must
_type be either "plain", "block" or "stationary" "plain"
block Block length for the block bootstrap,
_length mean block length for the stationary
bootstrap. Calculated by np::b.star()
if not provided NULL
n_rep Number of bootstrap replicates 1000

The only difference when testing for serial independence within a time series {Xt} is to create a two-column data set consisting of Xt and Xtk. For example, if we wish to perform this test for k=1 for one of the variables in the stock-exchange series, create the matrix of observations as below, and proceed exactly as in the i.i.d. case.

R> returns2 <- stock_data %>% select(SP500) %>% 
R+   mutate(sp500_lagged = lag(SP500))

Finally, the only thing that we must alter in order to perform the third test for serial cross-dependence is the bootstrap method. In the applications above, it suffices to use the standard bootstrap, where we resample with replacement from the data. This is implemented in the ind_test()-function by setting the bootstrap_type-argument to "plain", which is the default option. When testing for serial cross-dependence, we need to use a block-bootstrap procedure, and consider two options here: The block bootstrap with either fixed or random block sizes. This choice is specified by choosing bootstrap_type = "block" or bootstrap_type = "stationary", respectively, in the call to ind_test(). do not report significant differences in test performance using the different bootstrap types.

Test for financial contagion

Assume that we observe two financial time series {Xt} and {Yt} at times t=1,,T, and that some crisis occurs at time T<T. As described in Section 2, measure the local correlation between {Xt} and {Yt} before and after T, and take significant differences between these measurements as evidence against the null hypothesis of no linkage, or contagion, between the time series. In order to implement this test using the lg package, one must create two lg-objects: one for the observations covering the non-crisis period and one covering the crisis period. We do not enter a discussion here how to empirically identify such time periods; this is a job that must be done by the practitioner before performing the statistical test.

Let us illustrate the implementation of this test by looking at the same financial returns data that we have used in preceding sections. However, this time we will, in the spirit of , concentrate on GARCH(1,1)-filtrated daily returns on the S&P500 and FTSE100 indices from 2 January 1985 to 29 April 1987 in order to test for financial contagion between the US and UK stock markets following the global stock market crash of 19 October 1987 (“Black Monday”). Assume that these observations are loaded into the R workspace as the n1×2 data frame x_nc containing the observations covering the n1=728 days preceding the crisis and the n2×2 data frame x_c containing the observations covering the n2=140 consecutive trading days starting on Black Monday (see the online code supplement for details concerning the GARCH-filtration and data processing). In the code below, we construct one lg-object for each of these data frames with configuration matching the setup used by and perform the test by means of the cont_test()-function.

This function returns a list containing the estimated p-value as well as other useful statistics, including the empirical local correlations measured in the two time periods. See Table 4 for details concerning other arguments that may be passed to this function.

Table 4: Arguments to the cont_test()-function
Argument Explanation Default value
lg_object_nc The lg-object covering the non-crisis period
lg_object_c The lg-object covering the crisis period
grid_range Range of diagonal for measuring the LGC (5%, 95%) quantiles
grid_length The number of grid points to use 30
n_rep Number of bootstrap replicates 1000
weight Weight function function(y)
rep(1, nrow(y))
R> # Create the two lg-objects
R> lg_object_nc <- lg_main(x_nc, 
R+                         est_method = "5par",
R+                         transform_to_marginal_normality = FALSE)
R> 
R> lg_object_c <- lg_main(x_c, 
R+                         est_method = "5par",
R+                         transform_to_marginal_normality = FALSE)
R> 
R> # Run the test with a limited number of bootstrap replicates for 
R> # demonstration purposes.
R> result <- cont_test(lg_object_nc, lg_object_c, n_rep = 100)
R> 
R> # Print out the p-value
R> result$p_value

[1] 0.01

The small p-value means that we reject the null-hypothesis of no financial contagion between the time series after the crisis.

Partial local correlation

Consider the work finally by , who take the off-diagonal element in the local correlation matrix corresponding to the local conditional covariance matrix (12) or (14) as a local measure of conditional dependence between two stochastic variables X(1) and X(2) given the stochastic vector X(3). Furthermore, in the case with data having been transformed to marginal standard normality, they take the statistic (17)Tn,bCI=h(α^(z))dFn(z) as a measure of global conditional dependence. The lg package provides two key functions in this framework. The first, partial_cor(), calculates the local partial correlations as well as their estimated standard deviations on a specified grid in the (x(1),x(2),x(3))-space, and is essentially a wrapper function for the clg()-function presented in Section 4.1. See Table 5 for details. The second function, ci_test(), performs a test for conditional independence between the first two variables in a data set given the remaining variables using the test statistic ((17)) and a special bootstrap procedure (described briefly below) for approximating the null distribution.

Table 5: Arguments to the partial_cor()-function
Argument Explanation Default value
lg_object The lg-object created by lg_main()
grid The evaluation points for the LGPC, must be a
data frame or matrix having 2 columns NULL
cond Vector with fixed values for X(3) NULL
level Significance level for approximated confidence bands 0.95
Table 6: Arguments to the ci_test()-function
Argument Explanation Default value
lg_object The lg-object created by lg_main()
h The function h() in ((16)) function(x) x``^``2
n_rep Number of bootstrap replicates 500

It is well known in the econometrics literature that conditional independence tests are instrumental in the empirical detection of Granger causality . For example, if we continue to concentrate on the monthly stock returns data that we have already loaded into memory, we may test whether (18)H0:RFTSE100,tRSP500,t1|RFTSE100,t1 in the period starting in January 2009, the converse of which is a sufficient, but not necessary, condition for RSP500,t Granger causing RFTSE100,t. We perform the test by running the code below, where x is a data frame having the following columns strictly ordered as RFTSE100,t, RSP500,t1, and RFTSE100,t1 (see the online supplement for the pre-processing of data).

The critical values of this test are calculated using the bootstrap under the null hypothesis by independently resampling replicates from the conditional density estimates f^X(1)|X(3)(x(1)|x(3)) and f^X(2)|X(3)(x(2)|x(3)), as obtained by the clg()-function, using an approximated accept-reject algorithm. In order to avoid excessive optimization of the local likelihood function (1), we estimate fX(1)|X(3) and fX(2)|X(3) on the univariate regular grids x0(1) and x0(2), respectively (while keeping x(3) fixed at the observed values of X(3)), and produce interpolating functions f~X(1)|X(3) and f~X(2)|X(3) using cubic splines. It is much less computationally intensive to generate replicates from f~ than directly from f^.

We refer to the documentation of the lg package for details on how to finely tune the behavior of the bootstrapping algorithm by altering the arguments of the ci_test()-function and limiting our treatment to describing the arguments most suitable for modifications by the user in Table 6.

R> # Create the lg-object
R> lg_object <- lg_main(returns4)
R> 
R> # Perform the test
R> test_result <- ci_test(lg_object, n_rep = 100)
R> 
R> # Print out result
R> test_result$p_value

[1] 0.51

The conditional independence test does not provide evidence against the null-hypothesis (18).

5 Graphics

We conclude this article by describing the corplot() function for drawing dependence maps such as those displayed in Figure 1. report on such capabilities in the localgauss package, but the possibility of creating dependence maps was unfortunately removed from localgauss in the latest version 0.4.0 due to incompatibilities with the ggplot2 plotting engine. We make up for this loss by providing corplot(), a function that plots the estimated local correlations as provided by dlg(), or the estimated local partial correlations as provided by partial_cor().

The plotting function is highly customizable and provides a number of options covering most basic graphical options. Users well versed in the ggplot2 package may also modify the graphical object returned by corplot() in the standard way by adding layers as demonstrated in the example below.

In the first example, we generate a set of bivariate normally distributed data using the mvtnorm package and estimate the local Gaussian correlation on a regular grid using the dlg()-function. Passing the resulting dlg_object to corplot() without further arguments results in Figure 4.

R> # Make a regular grid in the domain of the distribution
R> grid <- expand.grid(seq(-3, 3, length.out = 7),
R+                     seq(-3, 3, length.out = 7))
R> 
R> x <- mvtnorm::rmvnorm(500, sigma = matrix(c(1, rho, rho, 1), 2))
R> lg_object <- lg_main(x, 
R+                      est_method = "5par", 
R+                      transform_to_marginal_normality = FALSE,
R+                      plugin_constant_joint = 4)
R> dlg_object <- dlg(lg_object, grid = grid)
R> 
R> # Make a dependence map using default setup
R> corplot(dlg_object)

We may tweak the appearance of our dependence map by passing further arguments to corplot(). Some of the options are demonstrated in the code chunk below, in which we, for example, superimpose the observations (by setting plot_obs = TRUE) and preventing the estimated local correlations from being plotted in areas without data. The latter option is available through the argument plot_thres, which works by calculating a bivariate kernel density estimate f~(x1,x2) for the pair of variables in question and only allowing ρ^(x1,x2) to be plotted if f~(x1,x2)/maxf~()> plot_thres. Adding layers to a dependence map using the ordinary ggplot2 syntax works as well, which we demonstrate in Figure 5 by changing the ggplot2 theme.

The plotting function works in the same way when plotting the local partial correlations returned by partial_cor(), and the arguments of corplot() are summarized in Table 7.

Table 7: Arguments to the corplot()-function
Argument Explanation Default value
dlg_object The object created by dlg() or partial_cor()
pair Which pair to plot if more than two variables 1
gaussian_scale Logical. Plot on the marginal st. normal scale? FALSE
plot_colormap Logical. Plot the colormap? TRUE
plot_obs Logical. Superimpose observations? FALSE
plot_labels Logical. Plot labels on dependence map? TRUE
plot_legend Logical. Add legend? FALSE
plot_thres Threshold for plotting the estimated LGC 0
alpha_tile Transparency of color tiles 0.8
alpha_point Transparency of points 0.8
low_color Color representing ρ^=1 "blue"
high_color Color representing ρ^=+1 "red"
break_int Break interval for color coding 0.2
label_size Size of labels in plot 3
font_family Font family for labels "sans"
point_size Size of points, if plotted NULL
xlim, ylim Axis limits NULL
xlab, ylab Axis labels NULL
rholab Title of legend NULL
main, subtitle Title and subtitle of plot NULL
R> corplot(dlg_object1, 
R+         plot_obs = TRUE, 
R+         plot_thres = 0.01,
R+         plot_labels = FALSE,
R+         alpha_point = 0.3,
R+         main = "",
R+         xlab = "",
R+         ylab = "") +
R+   theme_classic()
graphic without alt text
Figure 4: Dependence map produced by corplot() using the default configuration
graphic without alt text
Figure 5: Dependence map produced by tweaking the arguments of corplot()

6 Conclusion

The statistical literature has seen a number of applications of local Gaussian approximations in the last decade, covering several topics in dependence modeling and inference, as well as the estimation of multivariate density and conditional density functions. In this paper, we demonstrate the implementation of these methods in the R programming language using the lg package, as well as the graphical representation of the estimated local Gaussian correlation. The package is complete in the sense that all major methods that have been published within this framework is now easily accessible to the practitioner. The package is also designed with a modular infrastructure that allows future methodological developments using local Gaussian approximations to be easily added to the package.

7 Acknowledgements

The author gratefully acknowledges the constructive comments made by two anonymous referees.

CRAN packages used

lg, localgauss, magrittr, dplyr, ggplot2, mvtnorm

CRAN Task Views implied by cited packages

Databases, Distributions, Finance, ModelDeployment, Phylogenetics, Spatial, TeachingStatistics

Note

This article is converted from a Legacy LaTeX article using the texor package. The pdf version is the official version. To report a problem with the html, refer to CONTRIBUTE on the R Journal homepage.

Footnotes

    References

    S. M. Bache and H. Wickham. Magrittr: A forward-pipe operator for R. 2014. URL https://CRAN.R-project.org/package=magrittr. R package version 1.5.
    G. D. Berentsen, T. Kleppe and D. Tjøstheim. Introducing localgauss, an R-package for estimating and visualizing local Gaussian correlation. Journal of Statistical Software, 56(12): 1–18, 2014. URL https://doi.org/10.18637/jss.v056.i12.
    G. D. Berentsen and D. Tjøstheim. Recognizing and visualizing departures from independence in bivariate data using local Gaussian correlation. Statistics and Computing, 24(5): 785–801, 2014. URL https://doi.org/10.1007/s11222-013-9402-8.
    Datastream. 2018.
    K. J. Forbes and R. Rigobon. No contagion, only interdependence: Measuring stock market comovements. The Journal of Finance, 57(5): 2223–2261, 2002. URL https://doi.org/10.1111/0022-1082.00494.
    A. Genz, F. Bretz, T. Miwa, X. Mi, F. Leisch, F. Scheipl and T. Hothorn. Mvtnorm: Multivariate normal and t distributions. 2018. URL https://CRAN.R-project.org/package=mvtnorm. R package version 1.0-8.
    C. W. Granger. Testing for causality: A personal viewpoint. Journal of Economic Dynamics and control, 2: 329–352, 1980. URL https://doi.org/10.1016/0165-1889(80)90069-X.
    T. Hayfield and J. S. Racine. Nonparametric econometrics: The np package. Journal of Statistical Software, 27(5): 2008. URL https://doi.org/10.3929/ethz-b-000073064.
    N. L. Hjort and M. C. Jones. Locally parametric nonparametric density estimation. Annals of Statistics, 24(4): 1619–1647, 1996. URL https://doi.org/10.1214/aos/1032298288.
    R. A. Johnson and D. W. Wichern. Applied multivariate statistical analysis, sixth edition. Pearson Education Iternational, 2007.
    L. A. Jordanger. localgaussSpec: An r-package for local Gaussian spectral analysis. 2018. URL https://github.com/LAJordanger/localgaussSpec. R package.
    L. A. Jordanger and D. Tjøstheim. Nonlinear spectral analysis: A local gaussian approach. Journal of the American Statistical Association, 1–55, 2020.
    H. R. Kunsch. The jackknife and the bootstrap for general stationary observations. The Annals of Statistics, 17(3): 1217–1241, 1989. URL https://doi.org/10.1214/aos/1176347265.
    V. Lacal and D. Tjøstheim. Estimating and testing nonlinear local dependence between two time series. Journal of Business & Economic Statistics, 1–13, 2018.
    V. Lacal and D. Tjøstheim. Local Gaussian autocorrelation and tests of serial dependence. Journal of Time Series Analysis, 38(1): 51–71, 2017.
    H. Otneim. Lg: Locally gaussian distributions: Estimation and methods. 2019. URL https://CRAN.R-project.org/package=lg. R package version 0.4.0.
    H. Otneim. Multivariate and conditional density estimation using local Gaussian approximations. 2016.
    H. Otneim and D. Tjøstheim. Conditional density estimation using the local Gaussian correlation. Statistics and Computing, 28(2): 303–321, 2018.
    H. Otneim and D. Tjøstheim. The locally Gaussian density estimator for multivariate data. Statistics and Computing, 27(6): 1595–1616, 2017.
    H. Otneim and D. Tjøstheim. The locally Gaussian partial correlation. Journal of Business & Economic Statistics, 1–33, 2021.
    D. N. Politis and J. P. Romano. The stationary bootstrap. Journal of the American Statistical Association, 89(428): 1303–1313, 1994.
    B. W. Silverman. Density estimation for statistics and data analysis. Chapman; Hall, London, 1986.
    B. Støve, D. Tjøstheim and K. Hufthammer. Using local Gaussian correlation in a nonlinear re-examination of financial contagion. Journal of Empirical Finance, 25: 785–801, 2014.
    D. Tjøstheim and K. O. Hufthammer. Local Gaussian correlation: A new measure of dependence. Journal of Econometrics, 172(1): 33–48, 2013.
    H. Wickham. ggplot2: Elegant graphics for data analysis. Springer-Verlag New York, 2016. URL http://ggplot2.org.
    H. Wickham, R. François, L. Henry and K. Müller. Dplyr: A grammar of data manipulation. 2018. URL https://CRAN.R-project.org/package=dplyr. R package version 0.7.6.

    Reuse

    Text and figures are licensed under Creative Commons Attribution CC BY 4.0. The figures that have been reused from other sources don't fall under this license and can be recognized by a note in their caption: "Figure from ...".

    Citation

    For attribution, please cite this work as

    Otneim, "lg: An R package for Local Gaussian Approximations", The R Journal, 2021

    BibTeX citation

    @article{RJ-2021-079,
      author = {Otneim, Håkon},
      title = {lg: An R package for Local Gaussian Approximations},
      journal = {The R Journal},
      year = {2021},
      note = {https://rjournal.github.io/},
      volume = {13},
      issue = {2},
      issn = {2073-4859},
      pages = {38-56}
    }