# A Single or Multi-Parameter Optimization Engine (Beta Version) 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. utf8 ```r copBasic.fitpara.beta(uv=NULL, popstat=NULL, statf=NULL, cop=NULL, paradim=1, interval=NULL, par.init=NULL, ...) ``` ## Arguments - `uv`: An 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 or otherwise the `optim()` function in 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. ## Returns 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 **--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 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 to transform the initial guess $\alpha = \beta = 1/2$ into a domain more easily handled by the optimizer (`optim()` function in ). 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 ```r # See the Note section ```

copBasic.fitpara function