blomatrixCOP: A Matrix of Blomqvist-like Betas of a Copula
In wasquith/copBasic: General Bivariate Copula Theory and Many Utility Functions
blomatrixCOP R Documentation
A Matrix of Blomqvist-like Betas of a Copula
Description
Compute the Blomqvist-like Betas matrix \beta^\circ_\mathbf{C}-matrix of a copula, which is defined at presumably strategic points within \mathcal{I}^2, as (for as.blomCOPss=FALSE argument)
\beta^\circ_\mathbf{C} = \frac{\mathbf{C}(u^\circ,v^\circ)}{\mathbf{\Pi}(1/2, 1/2)} - 1\mbox{,}
where the u^\circ and v^\circ are of two types of gridded locations in \mathcal{I}^2 space and if u^\circ = 1/2 and v^\circ = 1/2, then central location of the matrix is Blomqvist Beta (blomCOP). The definition of \beta^\circ_\mathbf{C} is such that the matrix is entirely zero for the independence copula (\mathbf{\Pi}(u,v)) (P) when \mathbf{C}(u^\circ,v^\circ) = \mathbf{\Pi}(u,v) at the medial location u,v=1/2. Also, the definition here might be unique to the copBasic package. The decile version (blomatrixCOPdec) of this function uses u^\circ \in (1, 5, 9) / 10 and v^\circ \in (1, 5, 9) / 10. Whereas, the quartile version (blomatrixCOPiqr) of this function uses u^\circ \in (25, 50, 75) / 100 and v^\circ \in (25, 50, 75)/100. If as.blomCOPss=TRUE argument is set (default operation), then the coordinate locations in the matrix become the \beta^\diamond_\mathbf{C} of blomCOPss. As a rule \beta^\circ_\mathbf{C} \ne \beta^\diamond_\mathbf{C}.
Usage
blomatrixCOPdec(cop=NULL, para=NULL, as.sample=FALSE, as.blomCOPss=TRUE,
ctype=c("weibull", "hazen", "1/n",
"bernstein", "checkerboard"), ...)
blomatrixCOPiqr(cop=NULL, para=NULL, as.sample=FALSE, as.blomCOPss=TRUE,
ctype=c("weibull", "hazen", "1/n",
"bernstein", "checkerboard"), ...)
Arguments
cop
A copula function;
para
Vector of parameters or other data structure, if needed, to pass to the copula;
as.sample
A logical controlling whether an optional R data.frame in para is used to compute the \hat\beta^\circ_\mathbf{C}-matrix at which point the ctype argument will be passed to multiple calls of EMPIRcop;
as.blomCOPss
A logical to trigger blomCOPss for each of the (u,v) locations where for the blomCOPss calls (\beta^\diamond_\mathbf{C}(\mathbf{u}, \mathbf{v})): u \mapsto (u,u) \mapsto \mathbf{u} and v \mapsto (v,v) \mapsto \mathbf{v};
ctype
Argument of the same as EMPIRcop; and
...
Additional arguments to pass to the copula.
Value
The matrix for \beta^\circ_\mathbf{C} is returned depending on whether the decile or quartile version has been called.
Author(s)
W.H. Asquith
See Also
blomCOP, blomCOPss
Examples
round(blomatrixCOPdec(cop=P), digits=8); round(blomatrixCOPiqr(cop=P), digits=8)
# U|V=0.10 U|V=0.50 U|V=0.90 # U|V=0.25 U|V=0.50 U|V=0.75
# U|V=0.90 0 0 0 # U|V=0.75 0 0 0
# U|V=0.50 0 0 0 # U|V=0.50 0 0 0
# U|V=0.10 0 0 0 # U|V=0.25 0 0 0
round(blomatrixCOPdec(cop=PSP, as.blomCOPss=TRUE), digits=8)
# U|V=0.10 U|V=0.50 U|V=0.90
# U|V=0.90 0.4736842 0.8181818 0.5153268
# U|V=0.50 0.8181818 0.3333333 0.6459330
# U|V=0.10 0.8901099 0.4736842 0.1708292
round(blomatrixCOPdec(cop=PSP, as.blomCOPss=FALSE), digits=8)
# U|V=0.10 U|V=0.50 U|V=0.90
# U|V=0.90 0.09890110 0.81818182 4.2631579
# U|V=0.50 0.05263158 0.33333333 0.8181818
# U|V=0.10 0.01010101 0.05263158 0.0989011
## Not run:
set.seed(1)
td <- c(0.10, 0.50, 0.90)
UVn <- simCOP(n=5000, cop=glueCOP, col=8,
para=list(glue=0.4, cop1 =PLcop, cop2=PLcop,
para1=PLpar(rho=-0.5), para2=PLpar(rho=+0.5)))
points(td, rep(td[3], 3), cex=2, lwd=2, pch=3, col="red")
points(td, rep(td[2], 3), cex=2, lwd=2, pch=3, col="red")
points(td, rep(td[1], 3), cex=2, lwd=2, pch=3, col="red")
print(blomatrixCOPdec(as.sample=TRUE, para=UVn, ctype="weibull"))
# U|V=0.10 U|V=0.50 U|V=0.90
# U|V=0.90 -0.08222222 -0.580 -0.190
# U|V=0.50 0.30800000 0.112 0.264
# U|V=0.10 0.84000000 0.620 0.262
BMdn <- blomatrixCOPdec(cop=glueCOP,
para=list(glue=0.4, cop1=PLcop, cop2=PLcop,
para1=PLpar(rho=-0.5), para2=PLpar(rho=+0.5)))
print(round(BMdn, digits=8))
# U|V=0.10 U|V=0.50 U|V=0.90
# U|V=0.90 -0.08110464 -0.5569028 -0.2053772
# U|V=0.50 0.33449815 0.1202744 0.2424668
# U|V=0.10 0.75217766 0.6013719 0.2424668
## End(Not run)
## Not run:
set.seed(1); nsim <- 2000
para.pop <- list( cop1=GHcop, cop2=PLcop, alpha=0.359,
para1=c(4.003, 1.099), para2=0.882, beta=0.292)
UVs <- simCOP(nsim, cop=composite2COP, para=para.pop)
mtext("GIVEN THIS SAMPLE")
Rho <- rhoCOP(as.sample=TRUE, para=UVs) # Spearman Rho
BMn <- blomatrixCOPdec(as.sample=TRUE, para=UVs, ctype="weibull")
parafn <- function(k) {
c(exp(k[1])+1, exp(k[2]), exp(k[3]), pnorm(k[4]), pnorm(k[5]))
}
parafn_list <- function(k) {
k <- parafn(k)
list(cop1=GHcop, para1=c(k[1], k[2]), alpha=k[4],
cop2=PLcop, para2=k[4], beta=k[5])
}
BLOM_ofun <- function(para, statmat=NULL, parafn=NULL, rho=NA) {
para <- parafn(para)
new.para <- list(cop1=GHcop, para1=para[1:2], alpha=para[4],
cop2=PLcop, para2=para[3], beta=para[5])
bm <- blomatrixCOPdec(cop=composite2COP, para=new.para)
err <- sum((statmat - bm)^2) + (rhoCOP(cop=composite2COP, para=new.para) - rho)^2
#print(c(para, err))
return(err)
}
run1 <- function(graphics=TRUE, nsim=0) {
par.init <- c(log(1), log(1), log(1), qnorm(0.5), qnorm(0.5))
rt <- optim(par.init, BLOM_ofun, statmat=BMn, parafn=parafn, rho=Rho)
para <- parafn(rt$par)
para.fit <- list( cop1=GHcop, cop2=PLcop, alpha=para[4],
para1=para[1:2], para2=para[3], beta=para[5])
uv <- simCOP(nsim, cop=composite2COP, para=para.fit, graphics=graphics, col=2)
rmse <- round(rmseCOP(uv[,1], uv[,2], ctype="weibull",
cop=composite2COP, para=para.fit), digits=8)
if(graphics) mtext(paste0("RMSE(run1)=", rmse))
return(list(rmse=rmse, para=para.fit))
}
system.time(RUN1 <- run1(nsim=nsim))
par.init <- c(log(RUN1$para$para1[1]), log(RUN1$para$para1[2]),
log(RUN1$para$para2), qnorm(RUN1$para$alpha), qnorm(RUN1$para$beta))
RMSE_ofun <- function(para, parafn=NULL) {
para <- parafn(para)
new.para <- list(cop1=GHcop, para1=para[1:2], alpha=para[4],
cop2=PLcop, para2=para[3], beta=para[5])
new.rmse <- rmseCOP(UVs[,1], UVs[,2], cop=composite2COP, para=new.para)
#print(c(para, new.rmse))
return(new.rmse)
}
run2 <- function(graphics=TRUE, nsim=0, par.init=NULL) {
if(is.null(par.init)) {
par.init <- c(log(1), log(1), log(1), qnorm(0.5), qnorm(0.5))
}
rt <- optim(par.init, RMSE_ofun, parafn=parafn)
para <- parafn(rt$par)
para.fit <- list( cop1=GHcop, cop2=PLcop, alpha=para[4],
para1=para[1:2], para2=para[3], beta=para[5])
uv <- simCOP(nsim, cop=composite2COP, para=para.fit, graphics=graphics, col=4)
rmse <- round(rmseCOP(uv[,1], uv[,2], ctype="weibull",
cop=composite2COP, para=para.fit), digits=8)
if(graphics) mtext(paste0("RMSE(run2)=", rmse))
return(list(rmse=rmse, para=para.fit))
}
system.time(RUN2.1 <- run2(nsim=nsim, par.init=par.init))
system.time(RUN2.2 <- run2(nsim=nsim, par.init=NULL ))
GIVN <- c(para.pop$alpha, para.pop$para1, para.pop$beta, para.pop$para2)
FIT1 <- c(RUN1$para$alpha, RUN1$para$para1, RUN1$para$beta, RUN1$para$para2)
FIT1 <- round(FIT1, digits=3)
FIT2.1 <- c(RUN2.1$para$alpha, RUN2.1$para$para1, RUN2.1$para$beta, RUN2.1$para$para2)
FIT2.1 <- round(FIT2.1, digits=3)
FIT2.2 <- c(RUN2.2$para$alpha, RUN2.2$para$para1, RUN2.2$para$beta, RUN2.2$para$para2)
FIT2.2 <- round(FIT2.2, digits=3)
nms <- c("what", "alpha", "para1_1", "para1_2", "beta", "para2")
GIVN <- c("given", GIVN); FIT2.1 <- c("by_Blom1", FIT2.1)
FIT1 <- c("by_RMSE", FIT1); FIT2.2 <- c("by_Blom2", FIT2.2)
names(GIVN) <- nms; names(FIT1) <- nms
names(FIT2.1) <- nms; names(FIT2.2) <- nms
RESL <- cbind(data.frame(GIVN), data.frame(FIT1), data.frame(FIT2.1), data.frame(FIT2.2))
print(RESL) #
## End(Not run)
wasquith/copBasic documentation built on March 10, 2024, 11:24 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
| blomatrixCOP | R Documentation |
A Matrix of Blomqvist-like Betas of a Copula
Description
Compute the Blomqvist-like Betas matrix \beta^\circ_\mathbf{C}-matrix of a copula, which is defined at presumably strategic points within \mathcal{I}^2, as (for as.blomCOPss=FALSE argument)
\beta^\circ_\mathbf{C} = \frac{\mathbf{C}(u^\circ,v^\circ)}{\mathbf{\Pi}(1/2, 1/2)} - 1\mbox{,}
where the u^\circ and v^\circ are of two types of gridded locations in \mathcal{I}^2 space and if u^\circ = 1/2 and v^\circ = 1/2, then central location of the matrix is Blomqvist Beta (blomCOP). The definition of \beta^\circ_\mathbf{C} is such that the matrix is entirely zero for the independence copula (\mathbf{\Pi}(u,v)) (P) when \mathbf{C}(u^\circ,v^\circ) = \mathbf{\Pi}(u,v) at the medial location u,v=1/2. Also, the definition here might be unique to the copBasic package. The decile version (blomatrixCOPdec) of this function uses u^\circ \in (1, 5, 9) / 10 and v^\circ \in (1, 5, 9) / 10. Whereas, the quartile version (blomatrixCOPiqr) of this function uses u^\circ \in (25, 50, 75) / 100 and v^\circ \in (25, 50, 75)/100. If as.blomCOPss=TRUE argument is set (default operation), then the coordinate locations in the matrix become the \beta^\diamond_\mathbf{C} of blomCOPss. As a rule \beta^\circ_\mathbf{C} \ne \beta^\diamond_\mathbf{C}.
Usage
blomatrixCOPdec(cop=NULL, para=NULL, as.sample=FALSE, as.blomCOPss=TRUE,
ctype=c("weibull", "hazen", "1/n",
"bernstein", "checkerboard"), ...)
blomatrixCOPiqr(cop=NULL, para=NULL, as.sample=FALSE, as.blomCOPss=TRUE,
ctype=c("weibull", "hazen", "1/n",
"bernstein", "checkerboard"), ...)
Arguments
cop |
A copula function; |
para |
Vector of parameters or other data structure, if needed, to pass to the copula; |
as.sample |
A logical controlling whether an optional R |
as.blomCOPss |
A logical to trigger |
ctype |
Argument of the same as |
... |
Additional arguments to pass to the copula. |
Value
The matrix for \beta^\circ_\mathbf{C} is returned depending on whether the decile or quartile version has been called.
Author(s)
W.H. Asquith
See Also
blomCOP, blomCOPss
Examples
round(blomatrixCOPdec(cop=P), digits=8); round(blomatrixCOPiqr(cop=P), digits=8)
# U|V=0.10 U|V=0.50 U|V=0.90 # U|V=0.25 U|V=0.50 U|V=0.75
# U|V=0.90 0 0 0 # U|V=0.75 0 0 0
# U|V=0.50 0 0 0 # U|V=0.50 0 0 0
# U|V=0.10 0 0 0 # U|V=0.25 0 0 0
round(blomatrixCOPdec(cop=PSP, as.blomCOPss=TRUE), digits=8)
# U|V=0.10 U|V=0.50 U|V=0.90
# U|V=0.90 0.4736842 0.8181818 0.5153268
# U|V=0.50 0.8181818 0.3333333 0.6459330
# U|V=0.10 0.8901099 0.4736842 0.1708292
round(blomatrixCOPdec(cop=PSP, as.blomCOPss=FALSE), digits=8)
# U|V=0.10 U|V=0.50 U|V=0.90
# U|V=0.90 0.09890110 0.81818182 4.2631579
# U|V=0.50 0.05263158 0.33333333 0.8181818
# U|V=0.10 0.01010101 0.05263158 0.0989011
## Not run:
set.seed(1)
td <- c(0.10, 0.50, 0.90)
UVn <- simCOP(n=5000, cop=glueCOP, col=8,
para=list(glue=0.4, cop1 =PLcop, cop2=PLcop,
para1=PLpar(rho=-0.5), para2=PLpar(rho=+0.5)))
points(td, rep(td[3], 3), cex=2, lwd=2, pch=3, col="red")
points(td, rep(td[2], 3), cex=2, lwd=2, pch=3, col="red")
points(td, rep(td[1], 3), cex=2, lwd=2, pch=3, col="red")
print(blomatrixCOPdec(as.sample=TRUE, para=UVn, ctype="weibull"))
# U|V=0.10 U|V=0.50 U|V=0.90
# U|V=0.90 -0.08222222 -0.580 -0.190
# U|V=0.50 0.30800000 0.112 0.264
# U|V=0.10 0.84000000 0.620 0.262
BMdn <- blomatrixCOPdec(cop=glueCOP,
para=list(glue=0.4, cop1=PLcop, cop2=PLcop,
para1=PLpar(rho=-0.5), para2=PLpar(rho=+0.5)))
print(round(BMdn, digits=8))
# U|V=0.10 U|V=0.50 U|V=0.90
# U|V=0.90 -0.08110464 -0.5569028 -0.2053772
# U|V=0.50 0.33449815 0.1202744 0.2424668
# U|V=0.10 0.75217766 0.6013719 0.2424668
## End(Not run)
## Not run:
set.seed(1); nsim <- 2000
para.pop <- list( cop1=GHcop, cop2=PLcop, alpha=0.359,
para1=c(4.003, 1.099), para2=0.882, beta=0.292)
UVs <- simCOP(nsim, cop=composite2COP, para=para.pop)
mtext("GIVEN THIS SAMPLE")
Rho <- rhoCOP(as.sample=TRUE, para=UVs) # Spearman Rho
BMn <- blomatrixCOPdec(as.sample=TRUE, para=UVs, ctype="weibull")
parafn <- function(k) {
c(exp(k[1])+1, exp(k[2]), exp(k[3]), pnorm(k[4]), pnorm(k[5]))
}
parafn_list <- function(k) {
k <- parafn(k)
list(cop1=GHcop, para1=c(k[1], k[2]), alpha=k[4],
cop2=PLcop, para2=k[4], beta=k[5])
}
BLOM_ofun <- function(para, statmat=NULL, parafn=NULL, rho=NA) {
para <- parafn(para)
new.para <- list(cop1=GHcop, para1=para[1:2], alpha=para[4],
cop2=PLcop, para2=para[3], beta=para[5])
bm <- blomatrixCOPdec(cop=composite2COP, para=new.para)
err <- sum((statmat - bm)^2) + (rhoCOP(cop=composite2COP, para=new.para) - rho)^2
#print(c(para, err))
return(err)
}
run1 <- function(graphics=TRUE, nsim=0) {
par.init <- c(log(1), log(1), log(1), qnorm(0.5), qnorm(0.5))
rt <- optim(par.init, BLOM_ofun, statmat=BMn, parafn=parafn, rho=Rho)
para <- parafn(rt$par)
para.fit <- list( cop1=GHcop, cop2=PLcop, alpha=para[4],
para1=para[1:2], para2=para[3], beta=para[5])
uv <- simCOP(nsim, cop=composite2COP, para=para.fit, graphics=graphics, col=2)
rmse <- round(rmseCOP(uv[,1], uv[,2], ctype="weibull",
cop=composite2COP, para=para.fit), digits=8)
if(graphics) mtext(paste0("RMSE(run1)=", rmse))
return(list(rmse=rmse, para=para.fit))
}
system.time(RUN1 <- run1(nsim=nsim))
par.init <- c(log(RUN1$para$para1[1]), log(RUN1$para$para1[2]),
log(RUN1$para$para2), qnorm(RUN1$para$alpha), qnorm(RUN1$para$beta))
RMSE_ofun <- function(para, parafn=NULL) {
para <- parafn(para)
new.para <- list(cop1=GHcop, para1=para[1:2], alpha=para[4],
cop2=PLcop, para2=para[3], beta=para[5])
new.rmse <- rmseCOP(UVs[,1], UVs[,2], cop=composite2COP, para=new.para)
#print(c(para, new.rmse))
return(new.rmse)
}
run2 <- function(graphics=TRUE, nsim=0, par.init=NULL) {
if(is.null(par.init)) {
par.init <- c(log(1), log(1), log(1), qnorm(0.5), qnorm(0.5))
}
rt <- optim(par.init, RMSE_ofun, parafn=parafn)
para <- parafn(rt$par)
para.fit <- list( cop1=GHcop, cop2=PLcop, alpha=para[4],
para1=para[1:2], para2=para[3], beta=para[5])
uv <- simCOP(nsim, cop=composite2COP, para=para.fit, graphics=graphics, col=4)
rmse <- round(rmseCOP(uv[,1], uv[,2], ctype="weibull",
cop=composite2COP, para=para.fit), digits=8)
if(graphics) mtext(paste0("RMSE(run2)=", rmse))
return(list(rmse=rmse, para=para.fit))
}
system.time(RUN2.1 <- run2(nsim=nsim, par.init=par.init))
system.time(RUN2.2 <- run2(nsim=nsim, par.init=NULL ))
GIVN <- c(para.pop$alpha, para.pop$para1, para.pop$beta, para.pop$para2)
FIT1 <- c(RUN1$para$alpha, RUN1$para$para1, RUN1$para$beta, RUN1$para$para2)
FIT1 <- round(FIT1, digits=3)
FIT2.1 <- c(RUN2.1$para$alpha, RUN2.1$para$para1, RUN2.1$para$beta, RUN2.1$para$para2)
FIT2.1 <- round(FIT2.1, digits=3)
FIT2.2 <- c(RUN2.2$para$alpha, RUN2.2$para$para1, RUN2.2$para$beta, RUN2.2$para$para2)
FIT2.2 <- round(FIT2.2, digits=3)
nms <- c("what", "alpha", "para1_1", "para1_2", "beta", "para2")
GIVN <- c("given", GIVN); FIT2.1 <- c("by_Blom1", FIT2.1)
FIT1 <- c("by_RMSE", FIT1); FIT2.2 <- c("by_Blom2", FIT2.2)
names(GIVN) <- nms; names(FIT1) <- nms
names(FIT2.1) <- nms; names(FIT2.2) <- nms
RESL <- cbind(data.frame(GIVN), data.frame(FIT1), data.frame(FIT2.1), data.frame(FIT2.2))
print(RESL) #
## End(Not run)
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.