copBasic.fitpara: A Single or Multi-Parameter Optimization Engine (Beta...

In wasquith/copBasic: General Bivariate Copula Theory and Many Utility Functions

A Single or Multi-Parameter Optimization Engine (Beta Version)

Description

An example of a general implementation of a parameter optimization scheme using core features of the copBasic package. Because of the general complexity of the objectives for this function, the interface shown here is considered an “beta version” and nomenclature is subject to possibly sweeping changes in the future.

Usage

copBasic.fitpara.beta(uv=NULL, popstat=NULL, statf=NULL, cop=NULL,
                      paradim=1, interval=NULL, par.init=NULL, ...)

Arguments

uv	An R two column matrix or data.frame of a sample of nonexceedance probabilities u and v;
popstat	The population statistic(s) that will be used if uv is NULL;
statf	A function responsible at the minimum for computation of the theoretical values of the population statistic(s) that the optimization will revolve around; This function, if supporting an as.sample interface (e.g. hoefCOP) and if uv has been provided, will be dispatched to compute the population statistic(s);
cop	A copula function that is passed along to statf though of course the statf function can decide whether to use this argument or not;
paradim	The dimension of the parameters. In reality, the default triggers uni-dimensional root solving using the uniroot() function in R or otherwise the optim() function in R is used for multi-dimensional minimization with subtle changes in setup (see source code). Alternative logic could be have been used but it is felt that this is the most logical for future adaptations;
interval	The interval argument by the same name for the uniroot() function;
par.init	The initial parameter vector for the optim() function; and
...	Additional arguments to pass.

Value

A vector of the values for the parameter variable is returned

Note

One-Parameter Optimization—A demonstration of one-dimensional parameter optimimization using the Gini Gamma (giniCOP), which is a measure of bivariate association. There is no native support for Gini Gamma (and most of the other measures of association) in regards to being a parameter estimator at the copula function interface level in copBasic. (A converse example is one provided internally by the GHcop copula.)

  set.seed(24); n <- 230 # sample size to draw from Galambos copula but a
  # different copula was chosen for the fitting.
  sampleUV <- simCOP(n=n, cop=GLcop, para=1.5) # a random sample
  para <- copBasic.fitpara.beta(uv=sampleUV, statf=giniCOP,
                                interval=c(.1,200), cop=HRcop) # 1.959521

Three-Parameter Optimization—A demonstration of multi-dimensional parameter optimimization using the Gini Gamma (giniCOP), Nu-Skew (nuskewCOP), and Nu-Star (nustarCOP). This is substantially more complicated to implement. The Hüsler–Reiss copula (HRcop) is chosen both as part of the sample simulation for the study as well as the copula as part of the modeling. Using composition by the composite1COP, first establish a parent three parameter copula and simulate from it to create a bivariate sample in sampleUV that will be used for demonstration. A standard normal variate graphic of the simulation is generated by simCOP as well—later, additional results will be superimposed.

  n <- 610; set.seed(1999) # Seed value will be used again (see below)
  pop.para <- list(cop1=HRcop, para1=4, alpha=0.14, beta=0.35)
  sampleUV <- simCOP(n=n, cop=composite1COP, para=pop.para, col=3, snv=TRUE)

A custom objective function objfunc to serve as the statf for the copBasic.fitpara.beta call. The objective function has the as.sample interface (e.g. giniCOP) implemented for sample estimation. The most subtle feature of function presumably is the use of the pnorm() function in R for the \alpha and \beta parameters to keep each parameter in the range \alpha, \beta \in (0,1). Another subtly, which affects the choice of other copulas from HRcop, is how the parameter range of \Theta (the para[1] variable) is controlled—here the parameter remains untransformed but the lower domain edge is truncated by the return of infinity (return(Inf)). The getstat argument is only for short circuiting the objective so that objfunc can be used to compute the three statistics after the optimal parameters are found.

  "objfunc" <- function(para, stat=NA, as.sample=FALSE, getstat=FALSE, ...) {
      if(as.sample) {
         return(c(  giniCOP(para=para, as.sample=TRUE),
                  nuskewCOP(para=para, as.sample=TRUE),
                  nustarCOP(para=para, as.sample=TRUE)))
      }
      para[1]   <- ifelse(para[1] < 0, return(Inf), para[1]) # edge check
      para[2:3] <-  pnorm(para[2:3]) # detransform
      para <- list(cop1=HRcop, para1=para[1], alpha=para[2], beta=para[3])
      hp <- c(  giniCOP(composite1COP, para),
              nuskewCOP(composite1COP, para),
              nustarCOP(composite1COP, para))
      if(getstat) return(hp) # short circuit to get the statistics out
      dv <- stat; dv[dv == 0] <- 1 # changing dv "adapts" the error to
      return(sqrt(sum(((stat-hp)/dv)^2))) # trap division by zero
  }

The parameter estimation proceeds in the following code. The sample statistics (or target.stats) are computed and subsequently passed to the optimization using the popstat argument. Notice also the use of the qnorm() function in R to transform the initial guess \alpha = \beta = 1/2 into a domain more easily handled by the optimizer (optim() function in R). The transformed optimal parameters are computed, and then the values of the three statistics for the fit are computed. Lastly, a copBasic parameter object fit.para is created, which can be used for later comparisons.

  txt <- c("GiniGamma", "NuSkew", "NuStar")
  target.stats <- objfunc(sampleUV, as.sample=TRUE); names(target.stats) <- txt
  raw.fit.para <- copBasic.fitpara.beta(popstat=target.stats, statf=objfunc,
         par.init=c(1, qnorm(0.5), qnorm(0.5)), cop=composite1COP, paradim=3)
  fit.stats <- objfunc(raw.fit.para, getstat=TRUE); names(fit.stats) <- txt
  fit.para <- list(cop1=HRcop, para1=raw.fit.para[1],
                   alpha=pnorm(raw.fit.para[2]), beta=pnorm(raw.fit.para[3]))

The optimization is completed. It is informative to see what the simulation of the fitted copula looks like. Note: this particular example suffers from identifiability problems between the parameters—meaning that local minima exist or that satisfactory solutions using different parameters than used to generate the random sample can occur. The same seed is used so that one-to-one comparison of points can visually be made with the simCOP function call.

  set.seed(1999) # This value will be used again (see below)
  sampleUV <- simCOP(n=n, cop=composite1COP, para=fit.para,
                ploton=FALSE, pch=16, cex=0.5, col=2, snv=TRUE) # red dots

The visual comparison is qualitative. The tabular comparison of the sample statistics to those of the fitted model shows that perhaps an acceptable local minima has been attained in terms of “fit” but the subsequent comparison of the parameters of the population used to generate the sample and those of the fitted model seemingly depart substantially in the \alpha \rightarrow 0 parameter of the copula composition. The tail dependency parameters are similar, but further goodness-of-fit check is not made.

                       #   GiniGamma        NuSkew       NuStar
  print(target.stats)  #   0.5219027    -0.1940361    0.6108319
  print(fit.stats)     #   0.5182858    -0.1938848    0.6159566

                          # Parameter  Alpha       Beta
  print(ls.str(pop.para)) #      4.00   0.14      0.350  # given
  print(ls.str(fit.para)) #     11.2    0.187     0.427  # one solution

                                               # Tail Dependency Parameters
  taildepCOP(cop=composite1COP, para=pop.para) # lower=0 : upper=0.5838
  taildepCOP(cop=composite1COP, para=fit.para) # lower=0 : upper=0.5714(est.)

The demonstration ends with the comparison of the two asymmetrical density contours.

  densityCOPplot(cop=composite1COP, para=pop.para, contour.col=3)
  densityCOPplot(cop=composite1COP, para=fit.para, contour.col=2,
                                    ploton=FALSE)

Author(s)

W.H. Asquith

Examples

# See the Note section

wasquith/copBasic documentation built on March 10, 2024, 11:24 a.m.