In [6]:
library(VineCopula)
library(ggplot2)
In [7]:
dim <- 5
Matrix <- c(5, 2, 3, 1, 4,
0, 2, 3, 4, 1,
0, 0, 3, 4, 1,
0, 0, 0, 4, 1,
0, 0, 0, 0, 1)
Matrix <- matrix(Matrix, dim, dim)
# define R-vine pair-copula family matrix
family <- c(0, 1, 3, 4, 1,
0, 0, 3, 4, 1,
0, 0, 0, 4, 1,
0, 0, 0, 0, 3,
0, 0, 0, 0, 0)
family <- matrix(family, dim, dim)
# Params
par <- c(0, 0.2, 0.9, 1.5, 0.1,
0, 0, 1.1, 1.6, 0.3,
0, 0, 0, 1.9, 0.1,
0, 0, 0, 0, 4.8,
0, 0, 0, 0, 0)
par <- matrix(par, dim, dim)
# define second R-vine pair-copula parameter matrix
par2 <- matrix(0, dim, dim)
## define RVineMatrix object
RVM_0 <- VineCopula::RVineMatrix(Matrix = Matrix, family = family, par = par, par2 = par2)
RVM_0
In [147]:
n_theta <- 50
n <- 10000
medians <- rep(NA, n_theta)
means <- rep(NA, n_theta)
sums <- rep(NA, n_theta)
thetas <- seq(0.01, 0.99,length=n_theta)
par_new <- par
for (i in 1:n_theta){
par_new[5, 1] <- thetas[i]
RVM_new <- VineCopula::RVineMatrix(Matrix = Matrix, family = family, par = par_new, par2 = par2)
data <- VineCopula::RVineSim(n, RVM_new)
res <- VineCopula::RVineLogLik(data, RVM_0, separate = TRUE, calculate.V = FALSE)
medians[i] <- median(-res$loglik)
means[i] <- mean(-res$loglik)
sums[i] <- sum(-res$loglik)
}
In [148]:
df <- data.frame(thetas, medians, means, sums)
In [27]:
# Multiple plot function
#
# ggplot objects can be passed in ..., or to plotlist (as a list of ggplot objects)
# - cols: Number of columns in layout
# - layout: A matrix specifying the layout. If present, 'cols' is ignored.
#
# If the layout is something like matrix(c(1,2,3,3), nrow=2, byrow=TRUE),
# then plot 1 will go in the upper left, 2 will go in the upper right, and
# 3 will go all the way across the bottom.
#
multiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL) {
library(grid)
# Make a list from the ... arguments and plotlist
plots <- c(list(...), plotlist)
numPlots = length(plots)
# If layout is NULL, then use 'cols' to determine layout
if (is.null(layout)) {
# Make the panel
# ncol: Number of columns of plots
# nrow: Number of rows needed, calculated from # of cols
layout <- matrix(seq(1, cols * ceiling(numPlots/cols)),
ncol = cols, nrow = ceiling(numPlots/cols))
}
if (numPlots==1) {
print(plots[[1]])
} else {
# Set up the page
grid.newpage()
pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))
# Make each plot, in the correct location
for (i in 1:numPlots) {
# Get the i,j matrix positions of the regions that contain this subplot
matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))
print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row,
layout.pos.col = matchidx$col))
}
}
}
In [150]:
g_medians = ggplot(data=df, aes(x=thetas, y=medians)) +
geom_point(color="red")+ geom_smooth(color="grey40")+theme_bw()
g_means = ggplot(data=df, aes(x=thetas, y=means)) +
geom_point(color="red")+ geom_smooth(color="grey40")+theme_bw()
g_sums = ggplot(data=df, aes(x=thetas, y=sums)) +
geom_point(color="red")+ geom_smooth(color="grey40")+theme_bw()
multiplot(g_medians, g_means, g_sums, cols=1)
In [126]:
RVM_0
In [129]:
par_new[5, 1] <- thetas[means == min(means)]
RVM_new <- VineCopula::RVineMatrix(Matrix = Matrix, family = family, par = par_new, par2 = par2)
RVM_new
In [33]:
dim <- 2
Matrix <- c(2, 1,
0, 1)
Matrix <- matrix(Matrix, dim, dim)
# define R-vine pair-copula family matrix
family <- c(0, 3,
0, 0)
family <- matrix(family, dim, dim)
# Params
par <- c(0, 5.,
0, 0)
par <- matrix(par, dim, dim)
# define second R-vine pair-copula parameter matrix
par2 <- matrix(0, dim, dim)
## define RVineMatrix object
RVM_0 <- VineCopula::RVineMatrix(Matrix = Matrix, family = family, par = par, par2 = par2)
RVM_0
In [39]:
n_theta <- 20
n <- 10000
medians <- rep(NA, n_theta)
means <- rep(NA, n_theta)
sums <- rep(NA, n_theta)
grad <- rep(NA, n_theta)
thetas <- seq(1., 7,length=n_theta)
par_new <- par
for (i in 1:n_theta){
par_new[2, 1] <- thetas[i]
RVM_new <- VineCopula::RVineMatrix(Matrix = Matrix, family = family, par = par_new, par2 = par2)
data <- VineCopula::RVineSim(n, RVM_new)
res <- VineCopula::RVineLogLik(data, RVM_0, separate = TRUE, calculate.V = FALSE)
medians[i] <- median(res$loglik)
means[i] <- mean(res$loglik)
sums[i] <- sum(res$loglik)
grad[i] <- VineCopula::RVineGrad(data, RVM_0)$gradient**2
}
df <- data.frame(thetas, medians, means, sums)
In [41]:
g_medians = ggplot(data=df, aes(x=thetas, y=medians)) +
geom_point(color="red")+ geom_smooth(color="grey40")+theme_bw()
g_means = ggplot(data=df, aes(x=thetas, y=means)) +
geom_point(color="red")+ geom_smooth(color="grey40")+theme_bw()
g_sums = ggplot(data=df, aes(x=thetas, y=sums)) +
geom_point(color="red")+ geom_smooth(color="grey40")+theme_bw()
g_grad = ggplot(data=df, aes(x=thetas, y=grad)) +
geom_point(color="red")+ geom_smooth(color="grey40")+theme_bw()
multiplot(g_medians, g_means, g_sums, g_grad)
In [305]:
thetas[sums == min(sums)]
In [283]:
data = RVineSim(10000, RVM_0)
In [284]:
RVineMLE(data, RVM_new, maxit=1000)
In [11]:
sum(log(VineCopula::RVinePDF(data, RVM_target)))
In [23]:
library(VineCopula)
library(ggplot2)
dim <- 2
Matrix <- c(2, 1,
0, 1)
Matrix <- matrix(Matrix, dim, dim)
family <- c(0, 3,
0, 0)
family <- matrix(family, dim, dim)
par <- c(0, 5.,
0, 0)
par <- matrix(par, dim, dim)
par2 <- matrix(0, dim, dim)
RVM_target <- VineCopula::RVineMatrix(Matrix = Matrix, family = family, par = par, par2 = par2)
n_theta <- 20
n <- 10000
data_target <- VineCopula::RVineSim(n, RVM_target)
loglik <- rep(NA, n_theta)
grads <- rep(NA, n_theta)
thetas <- seq(0.01, 10.,length=n_theta)
par_new <- par
for (i in 1:n_theta){
par_new[2, 1] <- thetas[i]
RVM_new <- VineCopula::RVineMatrix(Matrix = Matrix, family = family, par = par_new, par2 = par2)
#data <- VineCopula::RVineSim(n, RVM_new)
res <- VineCopula::RVineLogLik(data_target, RVM_new, separate = FALSE, calculate.V = TRUE)
loglik[i] <- res$loglik
grads[i] <- VineCopula::RVineGrad(data, RVM_target)$gradient
}
df <- data.frame(thetas, loglik, grads)
ggplot(data=df, aes(x=thetas, y=loglik)) +
geom_point(color="red") +
geom_smooth(color="grey40") +
theme_bw()
ggplot(data=df, aes(x=thetas, y=grads**2)) +
geom_point(color="red") +
geom_smooth(color="grey40") +
theme_bw()
In [ ]:
In [ ]:
In [ ]: