Nothing
# Model: ADL-MIDAS (FADL-MIDAS) regression
# Estimation: MIDAS-NLS based on analytic derivatives (gradient and Hessian)
# Polynomial: Beta, 3 parametrizations, refer to "poly_spec"
# Paper, analytic derivatives:
# - Kostrov (2020) Estimating MIDAS regressions via MIDAS-NLS with revised optimization. Working paper.
# Optimizer: "nlminb", with constraints.
optim_ardl_beta <- function(y, z, x, poly_spec = 0, nbtrials = 100){
n <- length(y)
k <- dim(z)[2]
if (is.null(k)){
z <- matrix(z, nrow=length(z))
k <- 1
}
# append a vector of ones:
iota <- matrix(1, nrow = n, ncol = 1)
z <- cbind(iota, z)
# store data into a list:
dataList <- list (y=y, z=z, x=x, poly_spec=poly_spec)
# get parameter constraints:
constr <- get_constr_beta(poly_spec, k)
#-------------------- main optimization --------------------#
invisible(capture.output(opt <- multiStartoptim(objectivefn = estimate_ardl_beta,
gradient = gradient_ardl_beta,
data = dataList,
hessian = hessian_ardl_beta,
lower = constr$lw_b, upper = constr$up_b,
method = "nlminb",
nbtrials = nbtrials,
typerunif = "runifbase",
control = list(eval.max=1000, iter.max=1000,
rel.tol=1e-12, x.tol=1.5e-10))))
#-------------------- xxxxxxxxxxxxxxxx --------------------#
# back-out parameters:
k1 <- opt$par[1]
k2 <- opt$par[2]
k3 <- opt$par[3]
beta <- opt$par[4]
c <- opt$par[5]
rho <- opt$par[6:length(opt$par)]
# sort the output:
other_x <- NULL
for (i in 1:k)
other_x <- c(other_x, paste0("other-",i))
int <- "(Intercept)"
if (poly_spec == 0) {
coef <- matrix(c(c, rho, beta, k1, k2, k3), nrow = 1, ncol = k + 5)
rownames(coef) <- ""
colnames(coef) <- c(int,other_x,"beta","k1","k2", "k3")
} else if (poly_spec == 1) {
coef <- matrix(c(c, rho, beta, k2, k3), nrow = 1, ncol = k + 4)
rownames(coef) <- ""
colnames(coef) <- c(int,other_x,"beta","k2", "k3")
} else if (poly_spec == 2) {
coef <- matrix(c(c, rho, beta, k1, k2), nrow = 1, ncol = k + 4)
rownames(coef) <- ""
colnames(coef) <- c(int,other_x,"beta","k1", "k2")
}
class(coef) <- "midas.beta"
return(coef)
}
get_constr_beta <- function(poly_spec, k){
# define upper and lower constraints
up_theta1 <- 70
up_theta2 <- 70
up_theta3 <- 70
up_beta <- 40
up_rho <- 40
low_theta1 <- 0
low_theta2 <- 0
low_theta3 <- 0
low_beta <- -40
low_rho <- -40
# based on beta specification, construct vectors of constraints
if (poly_spec==0){
up_b = c(up_theta1, up_theta2, up_theta3, up_beta, rep(up_rho, k+1))
lw_b = c(low_theta1, low_theta2, low_theta3, low_beta, rep(low_rho, k+1))
} else if (poly_spec==1){
up_b = c(1, up_theta2, up_theta3, up_beta, rep(up_rho, k+1))
lw_b = c(1, low_theta2, low_theta3, low_beta, rep(low_rho, k+1))
} else if (poly_spec==2){
up_b = c(up_theta1, up_theta2, 0, up_beta, rep(up_rho, k+1))
lw_b = c(low_theta1, low_theta2, 0, low_beta, rep(low_rho, k+1))
} else if (poly_spec==3){
up_b = c(1, up_theta2, 0 , up_beta, rep(up_rho, k+1))
lw_b = c(1, low_theta2, 0, low_beta, rep(low_rho, k+1))
}
return(list(up_b = up_b, lw_b = lw_b))
}
estimate_ardl_beta <- function(args, dataList) {
y = dataList$y
z = dataList$z
x = dataList$x
p = dim(x)[2]
ii = matrix(1, p,1)
xx = matrix(c(1:p)/(p+1))
poly_spec = dataList$poly_spec
if (poly_spec==0){
theta1 = args[1]
theta3 = args[3]
} else if (poly_spec==1){
theta1 = 1
theta3 = args[3]
} else if (poly_spec==2){
theta1 = args[1]
theta3 = 0
} else if (poly_spec==3){
theta1 = 1
theta3 = 0
}
theta2 = args[2]
beta = args[4]
rho = args[5:length(args)]
m = gamma(theta1+theta2)/(gamma(theta1)*gamma(theta2))
weights = xx^(theta1-1)*(ii-xx)^(theta2-1)*m + theta3
#browser()
mse = t(y - z %*% rho - beta * x %*% weights)%*%(y - z %*% rho - beta * x %*% weights)
return(as.numeric(mse))
}
gradient_ardl_beta <- function(args, dataList) {
y = dataList$y
z = dataList$z
x = dataList$x
p = dim(x)[2]
ii = matrix(1, p,1)
xx = matrix(c(1:p)/(p+1))
poly_spec = dataList$poly_spec
if (poly_spec==0){
theta1 = args[1]
theta3 = args[3]
} else if (poly_spec==1){
theta1 = 1
theta3 = args[3]
} else if (poly_spec==2){
theta1 = args[1]
theta3 = 0
} else if (poly_spec==3){
theta1 = 1
theta3 = 0
}
theta2 = args[2]
beta = args[4]
rho = args[5:length(args)]
m = gamma(theta1+theta2)/(gamma(theta1)*gamma(theta2))
weights = xx^(theta1-1)*(ii-xx)^(theta2-1)*m + theta3
nabla_G = 2*(y- beta * x %*% weights - z%*%rho )
T_star_w_C = -beta * t(x)
T_star_beta_C = -t(weights) %*% t(x)
T_star_rho_C = - t(z)
weights_tilde = xx^(theta1-1)*(ii-xx)^(theta2-1)*m
log_derivative_theta2 = log(ii-xx) + digamma(theta1+theta2) - digamma(theta2)
T_star_theta2_W = t(log_derivative_theta2*weights_tilde)
part2 = T_star_theta2_W %*% T_star_w_C %*% nabla_G # theta2
if (poly_spec==0){
log_derivative_theta1 = log(xx) + digamma(theta1+theta2) - digamma(theta1)
T_star_theta1_W = t(log_derivative_theta1*weights_tilde)
T_star_theta3_W = t(ii)
part1 = T_star_theta1_W %*% T_star_w_C %*% nabla_G # theta1
part3 = T_star_theta3_W %*% T_star_w_C %*% nabla_G # theta3
} else if (poly_spec==1){
T_star_theta3_W = t(ii)
part1 = 0 # theta1
part3 = T_star_theta3_W %*% T_star_w_C %*% nabla_G # theta3
} else if (poly_spec==2){
log_derivative_theta1 = log(xx) + digamma(theta1+theta2) - digamma(theta1)
T_star_theta1_W = t(log_derivative_theta1*weights_tilde)
part1 = T_star_theta1_W %*% T_star_w_C %*% nabla_G # theta1
part3 = 0
} else if (poly_spec==3){
part1 = 0; # theta1
part3 = 0; # theta3
}
part4 = T_star_beta_C %*% nabla_G # beta
part5 = T_star_rho_C %*% nabla_G # rho
grad = c(part1, part2, part3, part4, part5)
return(grad)
}
hessian_ardl_beta <- function(args, dataList) {
y = dataList$y
z = dataList$z
x = dataList$x
p = dim(x)[2]
ii = matrix(1, p,1)
xx = matrix(c(1:p)/(p+1))
poly_spec = dataList$poly_spec
if (poly_spec==0){
theta1 = args[1]
theta3 = args[3]
} else if (poly_spec==1){
theta1 = 1
theta3 = args[3]
} else if (poly_spec==2){
theta1 = args[1]
theta3 = 0
} else if (poly_spec==3){
theta1 = 1
theta3 = 0
}
theta2 = args[2]
beta = args[4]
rho = args[5:length(args)]
k = matrix(0, 1, length(rho))
m = gamma(theta1+theta2)/(gamma(theta1)*gamma(theta2))
weights = xx^(theta1-1)*(ii-xx)^(theta2-1)*m + theta3
## line 1
E1 = 2*(y - beta * x %*% weights - z%*%rho )
F1 = -beta * t(x)
D1 = log(xx) + digamma(theta1+theta2) - digamma(theta1)
C1 = xx^(theta1-1)*(ii-xx)^(theta2-1)*m
### element 1
dC1_theta1 = C1 * D1
dD1_theta1 = trigamma(theta1+theta2 ) - trigamma(theta1)
dB1_theta1 = F1 %*% (-2*beta*x%*%dC1_theta1)
### element 2
D2 = log(ii-xx) + digamma(theta1+theta2) - digamma(theta2)
dC1_theta2 = C1*D2
dD1_theta2 = trigamma(theta1+theta2)
dB1_theta2 = F1%*%(-2*beta*x%*%dC1_theta2)
### element 3
dE1_theta3 = -2* beta*x %*% ii
### element 4
dF1_beta = -t(x)
dE1_beta = -2*x%*%weights
### element 5
dE1_rho = -2*z
## line 2
### element 2
dC2_theta2 = C1 * D2
dD2_theta2 = trigamma(theta1+theta2 ) - trigamma(theta2 )
### element 3
### element 4
### element 5
H22 = t(dC2_theta2 * D2 + C1 * dD2_theta2)%*%(F1%*%E1) + t(C1 *D2 )%*%dB1_theta2
H23 = t(C1 * D2)%*%F1%*%dE1_theta3
H24 = t(C1*D2) %*% (dF1_beta%*%E1+F1%*%dE1_beta)
H25 = t(C1 * D2)%*%F1%*%dE1_rho
## line 3
A3 = t(ii) %*% F1
### element 3
### element 4
dA3_beta = t(ii) %*% dF1_beta
### element 5
if (poly_spec==0){
H11 = t(dC1_theta1 * D1 + C1 * dD1_theta1)%*%(F1%*%E1) + t(C1 *D1 )%*%dB1_theta1
H12 = t(dC1_theta2 * D1 + C1 * dD1_theta2)%*%(F1%*%E1) + t(C1 *D1 )%*%dB1_theta2
H13 = t(C1 * D1)%*%F1%*%dE1_theta3
H14 = t(C1 *D1) %*% (dF1_beta%*%E1+F1%*%dE1_beta)
H15 = t(C1 * D1)%*%F1%*%dE1_rho
H33 = A3 %*% dE1_theta3
H34 = dA3_beta %*% E1 + A3 %*% dE1_beta
H35 = A3 %*% dE1_rho
} else if (poly_spec==1){
H11 = 0
H12 = 0
H13 = 0
H14 = 0
H15 = k
H33 = A3 %*% dE1_theta3
H34 = dA3_beta %*% E1 + A3 %*% dE1_beta
H35 = A3 %*% dE1_rho
} else if (poly_spec==2){
H11 = t(dC1_theta1 * D1 + C1 * dD1_theta1)%*%(F1%*%E1) + t(C1 *D1 )%*%dB1_theta1
H12 = t(dC1_theta2 * D1 + C1 * dD1_theta2)%*%(F1%*%E1) + t(C1 *D1 )%*%dB1_theta2
H13 = 0
H14 = t(C1 *D1) %*% (dF1_beta%*%E1+F1%*%dE1_beta)
H15 = t(C1 * D1)%*%F1%*%dE1_rho
H23 = 0
H33 = 0
H34 = 0
H35 = k
} else if (poly_spec==3){
H11 = 0
H12 = 0
H13 = 0
H14 = 0
H15 = k
H23 = 0
H33 = 0
H34 = 0
H35 = k
}
## line 4
A4 = -t(weights) %*% t(x)
### element 4
### element 5
H44 = A4 %*% dE1_beta
H45 = A4 %*% dE1_rho
## line 5
### element 5
H55 = -t(z) %*% dE1_rho
Hess = rbind( c(H11, H12, H13, H14, H15),
c(H12, H22, H23, H24, H25),
c(H13, H23, H33, H34, H35),
c(H14, H24, H34, H44, H45),
cbind(t(H15), t(H25), t(H35), t(H45), H55))
#browser()
return(Hess)
}
# Model: ADL-MIDAS (FADL-MIDAS) regression
# Estimation: MIDAS-NLS based on analytic derivatives (gradient and Hessian)
# Polynomial: Exponential Almon parametrized by theta1 and theta 2.
# Paper, analytic derivatives: based on
# - Kostrov A. and Tetereva A. (2018) Forecasting realized correlations: a MIDAS approach. Working paper.
# - https://www.researchgate.net/publication/331498346_Forecasting_realized_correlations_a_MIDAS_approach
# Optimizer: "nlminb" with constraints.
optim_ardl_expalmon <- function(y, z, x, nbtrials = 100){
n <- length(y)
k <- dim(z)[2]
if (is.null(k)){
z <- matrix(z, nrow=length(z))
k <- 1
}
# append a vector of ones:
iota <- matrix(1, nrow = n, ncol = 1)
z <- cbind(iota, z)
# store data into a list:
dataList <- list (y=y, z=z, x=x)
# get parameter constraints:
constr <- get_constr_expalmon(k)
#-------------------- main optimization --------------------#
invisible(capture.output(opt <- multiStartoptim(objectivefn = estimate_ardl_expalmon,
gradient = gradient_ardl_expalmon,
data = dataList,
hessian = hessian_ardl_expalmon,
lower = constr$lw_b, upper = constr$up_b,
method = "nlminb",
nbtrials = nbtrials,
typerunif = "runifbase",
control = list(eval.max=1000, iter.max=1000,
rel.tol=1e-12, x.tol=1.5e-10))))
#-------------------- xxxxxxxxxxxxxxxx --------------------#
# back-out parameters:
k1 <- opt$par[1]
k2 <- opt$par[2]
beta <- opt$par[3]
c <- opt$par[4]
rho <- opt$par[5:length(opt$par)]
# sort the output:
other_x <- NULL
for (i in 1:k)
other_x <- c(other_x, paste0("other-",i))
int <- "(Intercept)"
coef <- matrix(c(c, rho, beta, k1, k2), nrow = 1, ncol = k + 4)
rownames(coef) <- ""
colnames(coef) <- c(int,other_x,"beta","k1","k2")
class(coef) <- "midas.expalmon"
return(coef)
}
get_constr_expalmon <- function(k){
# define upper and lower constraints
up_beta <- 20
up_rho <- 20
up_theta1 <- 2
up_theta2 <- 0
low_beta <- -10
low_rho <- -20
low_theta1 <- 0
low_theta2 <- -2
up_b <- c(up_theta1 , up_theta2, up_beta, rep(up_rho, k+1))
lw_b <- c(low_theta1, low_theta2, low_beta, rep(low_rho, k+1))
return(list(up_b = up_b, lw_b = lw_b))
}
estimate_ardl_expalmon <- function(args, dataList) {
y <- dataList$y
z <- dataList$z
x <- dataList$x
p <- dim(x)[2]
xi <- matrix(1:p, p,1)
xi_sq <- xi^2
theta1 <- args[1]
theta2 <- args[2]
beta <- args[3]
rho <- args[4:length(args)]
weights <- exp(theta1 * xi + theta2 * xi_sq)
#browser()
mse <- t(y - z %*% rho - beta * x %*% weights)%*%(y - z %*% rho - beta * x %*% weights)
return(as.numeric(mse))
}
gradient_ardl_expalmon <- function(args, dataList) {
y <- dataList$y
z <- dataList$z
x <- dataList$x
p <- dim(x)[2]
xi <- matrix(1:p, p,1)
xi_sq <- xi^2
theta1 <- args[1]
theta2 <- args[2]
beta <- args[3]
rho <- args[4:length(args)]
weights <- exp(theta1 * xi + theta2 * xi_sq)
nabla_G <- 2*(y- beta * x %*% weights - z%*%rho )
T_star_w_C <- -beta * t(x)
T_star_beta_C <- -t(weights) %*% t(x)
T_star_rho_C <- - t(z)
T_star_theta1_W <- t(weights * xi)
T_star_theta2_W <- t(weights * xi_sq)
part1 <- T_star_theta1_W %*% T_star_w_C %*% nabla_G # theta1
part2 <- T_star_theta2_W %*% T_star_w_C %*% nabla_G # theta2
part3 <- T_star_beta_C %*% nabla_G # beta
part4 <- T_star_rho_C %*% nabla_G # rho
grad <- c(part1, part2, part3, part4)
return(grad)
}
hessian_ardl_expalmon <- function(args, dataList) {
y <- dataList$y
z <- dataList$z
x <- dataList$x
p <- dim(x)[2]
xi <- matrix(1:p, p,1)
xi_sq <- xi^2
theta1 <- args[1]
theta2 <- args[2]
beta <- args[3]
rho <- args[4:length(args)]
weights <- exp(theta1 * xi + theta2 * xi_sq)
## line 1
A1_1 <- t(weights * xi)
A1_2 <- -beta * t(x)
A1 <- A1_1 %*% A1_2
B1 <- 2*(y - beta * x %*% weights - z%*%rho )
### element 1
dA1_theta1 <- t(weights * xi_sq) %*% A1_2
dB1_theta1 <- - 2* beta * x %*% (weights * xi)
### element 2
dA1_theta2 <- t(weights * xi^3) %*% A1_2
dB1_theta2 <- - 2* beta * x %*% (weights * xi_sq)
### element 3
dA1_beta <- t(weights * xi)%*%t(-x)
dB1_beta <- -2*x%*%weights
### element 4
dB1_rho <- -2*z
H11 <- dA1_theta1 %*% B1 + A1 %*% dB1_theta1
H12 <- dA1_theta2 %*% B1 + A1 %*% dB1_theta2
H13 <- dA1_beta %*% B1 + A1 %*% dB1_beta
H14 <- A1 %*% dB1_rho
## line 2
A2_1 <- t(weights * xi_sq)
A2 <- A2_1 %*% A1_2
### element 2
dA2_theta2 <- t(weights * xi^4) %*% A1_2
### element 3
dA2_beta <- A2_1 %*% t(-x)
### element 4
H22 <- dA2_theta2 %*% B1 + A2 %*% dB1_theta2
H23 <- dA2_beta %*% B1 + A2 %*% dB1_beta
H24 <- A2 %*% dB1_rho
## line 3
A3 <- -t(weights) %*% t(x)
### element 3
### element 4
### element 5
H33 <- A3 %*% dB1_beta
H34 <- A3 %*% dB1_rho
## line 4
H44 <- t(-z) %*% dB1_rho
Hess <- rbind( c(H11, H12, H13, H14),
c(H12, H22, H23, H24),
c(H13, H23, H33, H34),
cbind(t(H14), t(H24), t(H34), H44))
#browser()
return(Hess)
}
# Model: MIDAS Logit regression
# Estimation: MIDAS-NLS based on analytic derivatives (gradient and Hessian)
# Polynomial: Exponential Almon parametrized by theta1 and theta 2
# Paper, analytic derivatives: based on
# - F. Audrino, A. Kostrov, and J.-P. Ortega (2019) Predicting U.S. Bank Failures with MIDAS Logit Models, JFQA, 54(6), 2575-2603.
# - https://doi.org/10.1017/S0022109018001308
# Optimizer: nlm.
optim_logit_expalmon <- function(y, z, x){
if (!all(unique(y)%in%c(0,1))){
stop("response variable must be binary for MIDAS Logit model.")
}
n <- length(y)
k <- dim(z)[2]
if (is.null(k)){
z <- matrix(z, nrow=length(z))
k <- 1
}
# append a vector of ones:
iota <- matrix(1, nrow = n, ncol = 1)
z <- cbind(iota, z)
# initial values
init = rep(0,3 + k+1)
#-------------------- main optimization --------------------#
opt <- suppressWarnings(stats::nlm(f = estimateLogit_expalm_nlm,
p = init,
y=y, V=z, Z=x,
hessian = TRUE,
check.analyticals = TRUE))
#-------------------- xxxxxxxxxxxxxxxx --------------------#
# back-out parameters:
k1 <- opt$estimate[1]
k2 <- opt$estimate[2]
beta <- opt$estimate[3]
c <- opt$estimate[4]
rho <- opt$estimate[5:length(opt$estimate)]
# sort the output:
fn <- NULL
for (i in 1:k)
fn <- c(fn, paste0("factor-",i))
int <- "(Intercept)"
coef <- matrix(c(c, rho, beta, k1, k2), nrow = 1, ncol = k + 4)
rownames(coef) <- ""
colnames(coef) <- c(int,fn,"beta","k1","k2")
class(coef) <- "logit.expalmon"
return(coef)
}
estimateLogit_expalm_nlm <- function(args, y, V, Z){
iota = matrix(1, length(y) )
p = dim(Z)[2]
xi = matrix(1:p, p,1)
xi_sq = xi^2
theta1 = args[1]
theta2 = args[2]
beta = args[3]
rho = args[4:length(args)]
weights = exp(theta1 * xi + theta2 * xi_sq)
LogLM = -( t(y) %*% (V %*% rho + beta * Z %*% weights) -
t(iota) %*% log(iota+exp(V %*% rho + beta * Z %*% weights)) )
out = as.numeric(LogLM)
if (is.na(out)==TRUE) {
out = Inf
}
attr(out, 'gradient') <- gradLogit_expalm(args, y, V, Z)
attr(out, 'hessian') <- HessianLogit_expalm(args, y, V, Z)
return(out)
}
gradLogit_expalm <- function(args, y, V, Z) {
p = dim(Z)[2]
xi = matrix(1:p, p,1)
xi_sq = xi^2
theta1 = args[1]
theta2 = args[2]
beta = args[3]
rho = args[4:length(args)]
weights = exp(theta1 * xi + theta2 * xi_sq)
nabla_G = - ( y- exp(beta * Z %*% weights + V%*%rho) *
1 / (1 + exp(beta * Z %*% weights + V%*%rho)))
T_star_w_C = beta * t(Z)
T_star_beta_C = t(weights) %*% t(Z)
T_star_rho_C = t(V)
T_star_theta1_W = t(weights * xi)
T_star_theta2_W = t(weights * xi_sq)
part1 = T_star_theta1_W %*% T_star_w_C %*% nabla_G # theta1
part2 = T_star_theta2_W %*% T_star_w_C %*% nabla_G # theta2
part3 = T_star_beta_C %*% nabla_G # beta
part4 = T_star_rho_C %*% nabla_G # rho
grad = c(part1, part2, part3, part4)
return(grad)
}
HessianLogit_expalm <- function(args, y, V, Z) {
p = dim(Z)[2]
xi = matrix(1:p, p,1)
xi_sq = xi^2
theta1 = args[1]
theta2 = args[2]
beta = args[3]
rho = args[4:length(args)]
weights = exp(theta1 * xi + theta2 * xi_sq)
lincomb = beta * Z %*% weights + V%*%rho
## line 1
A1_1 = t(weights * xi)
A1_2 = beta * t(Z)
A1 = A1_1 %*% A1_2
B1 = - ( y- exp(lincomb) * 1 / (1 + exp(lincomb)))
dB1_lincomb = exp(lincomb)*(1+exp(lincomb))^(-1) - exp(lincomb)^(2) * (1+exp(lincomb))^(-2)
### element 1
dA1_theta1 = t(weights * xi_sq) %*% A1_2
dB1_theta1 = dB1_lincomb * (beta * Z) %*% (weights * xi)
### element 2
dA1_theta2 = t(weights * xi^3) %*% A1_2
dB1_theta2 = dB1_lincomb * (beta * Z) %*% (weights * xi_sq)
### element 3
dA1_beta = t(weights * xi)%*%t(Z)
dB1_beta = dB1_lincomb * (Z%*%weights)
### element 4
dB1_rho = diag( as.vector(dB1_lincomb) ) %*% V
H11 = dA1_theta1 %*% B1 + A1 %*% dB1_theta1
H12 = dA1_theta2 %*% B1 + A1 %*% dB1_theta2
H13 = dA1_beta %*% B1 + A1 %*% dB1_beta
H14 = A1 %*% dB1_rho
## line 2
A2_1 = t(weights * xi_sq)
A2 = A2_1 %*% A1_2
### element 2
dA2_theta2 = t(weights * xi^4) %*% A1_2
### element 3
dA2_beta = A2_1 %*% t(Z)
### element 4
H22 = dA2_theta2 %*% B1 + A2 %*% dB1_theta2
H23 = dA2_beta %*% B1 + A2 %*% dB1_beta
H24 = A2 %*% dB1_rho
## line 3
A3 = t(weights) %*% t(Z)
### element 3
### element 4
### element 5
H33 = A3 %*% dB1_beta
H34 = A3 %*% dB1_rho
## line 4
H44 = t(V) %*% dB1_rho
Hess = rbind( c(H11, H12, H13, H14),
c(H12, H22, H23, H24),
c(H13, H23, H33, H34),
cbind(t(H14), t(H24), t(H34), H44))
#browser()
return(Hess)
}
optim_ardl_legendre <- function(y, z, x, legendre_degree = 3){
n <- length(y)
k <- dim(z)[2]
if (is.null(k)){
z <- matrix(z, nrow=length(z))
k <- 1
}
# append a vector of ones:
iota <- matrix(1, nrow = n, ncol = 1)
z <- cbind(iota, z)
p <- lb(degree = legendre_degree, jmax = dim(x)[2])
x <- x%*%p
bigx <- cbind(z, x)
coef <- t(solve(t(bigx)%*%bigx)%*%t(bigx)%*%y)
other_x <- NULL
for (i in 1:k)
other_x <- c(other_x, paste0("other-",i))
poly_n <- NULL
for (i in 0:legendre_degree)
poly_n <- c(poly_n, paste0("Poly-",i))
int <- "(Intercept)"
rownames(coef) <- ""
colnames(coef) <- c(int,other_x,poly_n)
class(coef) <- "midas.legendre"
return(coef)
}
optim_logit_legendre <- function(y, z, x, legendre_degree = 3){
n <- length(y)
k <- dim(z)[2]
if (is.null(k)){
z <- matrix(z, nrow=length(z))
k <- 1
}
# append a vector of ones:
iota <- matrix(1, nrow = n, ncol = 1)
z <- cbind(iota, z)
p <- lb(degree = legendre_degree, jmax = dim(x)[2])
x <- x%*%p
bigx <- cbind(z, x)
coef <- stats::glm(y~-1+bigx)$coef
coef <- matrix(coef, ncol = length(coef))
fn <- NULL
for (i in 1:k)
fn <- c(fn, paste0("factor-",i))
poly_n <- NULL
for (i in 0:legendre_degree)
poly_n <- c(poly_n, paste0("Poly-",i))
int <- "(Intercept)"
rownames(coef) <- ""
colnames(coef) <- c(int,fn,poly_n)
class(coef) <- "logit.legendre"
return(coef)
}
multiStartoptim <- function(start0 = NULL, objectivefn, gradient = NULL, ..., hessian = NULL,
lower = -Inf, upper = Inf, control = list(),
method = c("L-BFGS-B", "Nelder-Mead","nlminb"), nbtrials = NULL,
typerunif = c("runifbase", "runifantithetics", "sobol", "torus", "niederreiterlodisp"), localsearch = c("exhaustive", "median"),
verb=FALSE, nbclusters = NULL)
{
if(missing(objectivefn)) stop("Objective function must be provided")
nbpar <- length(lower)
if (nbpar != length(upper) || lower > upper) stop("parameter lower must be lower than parameter upper componentwise, and they must have the same length")
if (missing(method) || !(method %in% c("L-BFGS-B", "Nelder-Mead", "nlminb")))
{
warning("optimization method is missing or inadequate, default is set to nlminb")
themethod <- "nlminb"
}
else
{
themethod <- match.arg(method)
}
if (missing(nbclusters) || is.null(nbclusters) || nbclusters <= 1)
{
if (missing(start0) && (missing(nbtrials) || nbtrials <= 0)) stop("When nbtrials is missing or equal to 0, starting parameters 'start0' should be provided")
if (missing(nbtrials) || is.null(nbtrials) || nbtrials <= 0)
{
if (!missing(typerunif)) warning("unused argument typerunif")
if (missing(hessian) && themethod %in% c("L-BFGS-B", "Nelder-Mead")) hessian <- FALSE
return(call_localoptim(start0 = start0, objectivefn = objectivefn, gradient = gradient, ..., hessian = hessian,
lower = lower, upper = upper, control = control, themethod = themethod, nbtrials = nbtrials ,
verb = verb, cl = NULL))
}
else
{
if (!missing(start0)) warning("starting value is unused when nbtrials is provided")
if (missing(typerunif) || !(typerunif %in% c("runifbase", "runifantithetics", "sobol", "torus", "niederreiterlodisp")))
{
warning("typerunif is either missing or inadequate, and was set to default 'runifbase'")
typerunif <- "runifbase"
}
else
{
typerunif <- match.arg(typerunif)
}
if (!is.finite(lower) || !is.finite(upper))
stop("Both the lower and upper bounds for parameters should be provided and finite")
localsearch <- match.arg(localsearch)
if (length(localsearch) == 0) localsearch <- "exhaustive"
if (localsearch == "exhaustive")
{
start0 <- rUnif(nbtrials = nbtrials, nbpar = nbpar, lower = lower, upper = upper,
method=typerunif)
}
if (localsearch == "median")
{
if (nbpar == 1)
{
U <- rUnif(nbtrials = nbtrials, nbpar = nbpar, lower = lower, upper = upper,
method=typerunif)
fwU <- objectivefn(U, ...)
start0 <- U[fwU < median(fwU)]
nbtrials <- length(start0)
}
else
{
U <- rUnif(nbtrials = nbtrials, nbpar = nbpar, lower = lower, upper = upper,
method=typerunif)
fwU <- apply(U, 1, function(x) objectivefn(x, ...))
start0 <- U[fwU < median(fwU), ]
nbtrials <- dim(start0)[1]
}
}
if (missing(hessian) && themethod %in% c("L-BFGS-B", "Nelder-Mead")) hessian <- FALSE
return(call_localoptim(start0 = start0, objectivefn = objectivefn, gradient = gradient, ..., hessian = hessian,
lower = lower, upper = upper, control = control, themethod = themethod, nbtrials = nbtrials ,
verb = verb, cl = NULL))
}
}
else
{
if (nbclusters < 0 || floor(nbclusters) != nbclusters || is.array(nbclusters) || !is.numeric(nbclusters)) stop("nbclusters must be a positive integer")
if(missing(lower) || missing(upper)) stop("lower and upper are missing and must be provided")
if (!is.finite(lower) || !is.finite(upper)) stop("Both the lower and upper bounds for parameters should be provided and finite")
if(verb) warning("argument verb is unused in parallel computation")
if (missing(nbtrials) || is.null(nbtrials) || nbtrials <= 0)
{
warning("nbtrials is not provided (or negative), default is set to 50")
nbtrials <- 50
}
if (missing(typerunif) || !(typerunif %in% c("runifbase", "runifantithetics", "sobol", "torus", "niederreiterlodisp")))
{
warning("typerunif is either missing or inadequate, and was set to default 'runifbase'")
typerunif <- "runifbase"
} else
{
typerunif <- match.arg(typerunif)
}
localsearch <- match.arg(localsearch)
if (length(localsearch) == 0) localsearch <- "exhaustive"
if (localsearch == "exhaustive")
{
start0 <- rUnif(nbtrials = nbtrials, nbpar = nbpar, lower = lower, upper = upper,
method=typerunif)
}
if (localsearch == "median")
{
if (nbpar == 1)
{
U <- rUnif(nbtrials = nbtrials, nbpar = nbpar, lower = lower, upper = upper,
method=typerunif)
fwU <- objectivefn(U, ...)
start0 <- U[fwU < median(fwU)]
nbtrials <- length(start0)
}
else
{
U <- rUnif(nbtrials = nbtrials, nbpar = nbpar, lower = lower, upper = upper,
method=typerunif)
fwU <- apply(U, 1, function(x) objectivefn(x, ...))
start0 <- U[fwU < median(fwU), ]
nbtrials <- dim(start0)[1]
}
}
if (missing(hessian) && themethod %in% c("L-BFGS-B", "Nelder-Mead")) hessian <- FALSE
packageStartupMessage("Processing...", appendLF = FALSE)
cat("\n")
return(call_localoptim(start0 = start0, objectivefn = objectivefn, gradient = gradient, ..., hessian = hessian,
lower = lower, upper = upper, control = control, themethod = themethod, nbtrials = nbtrials ,
verb = verb, cl = nbclusters))
}
}
rUnif <- function(nbtrials, nbpar, lower, upper,
method=c("runifbase", "runifantithetics", "sobol", "torus", "niederreiterlodisp"))
{
callrunifmat <- function(nbtrials, nbpar, lower, upper)
{
if (nbpar == 1)
{
return(lower + (upper - lower) * randtoolbox::SFMT(nbtrials))
}
else
{
Umin <- matrix(rep.int(lower, nbtrials), nrow = nbtrials, ncol = nbpar, byrow=T)
etendue <- upper - lower
return(Umin + matrix(rep.int(etendue, nbtrials), nrow = nbtrials, ncol=nbpar, byrow=T)*randtoolbox::SFMT(nbtrials,nbpar))
}
}
callrunifmatantith <- function(nbtrials, nbpar, lower, upper)
{
nbtrials_2 <- 0.5*nbtrials
if (nbpar == 1)
{
Uni <- randtoolbox::SFMT(nbtrials_2)
return(lower + (upper - lower) * c(Uni, 1-Uni))
}
else
{
Umin <- matrix(rep.int(lower, nbtrials), nrow = nbtrials, ncol = nbpar, byrow=T)
etendue <- upper - lower
Uni <- randtoolbox::SFMT(nbtrials_2,nbpar)
return(Umin + matrix(rep.int(etendue, nbtrials), nrow = nbtrials, ncol = nbpar, byrow = T)*rbind(Uni, 1-Uni))
}
}
callxlodisp <- function(nbtrials, lower, upper)
{
nbpar <- length(lower)
if (nbpar == 1)
{
ylodisp <- log(2*seq_len(nbtrials)[-1] - 3)/(log(2))
xlodisp <- c(1, ylodisp - floor(ylodisp))
return(lower + (upper - lower) * xlodisp)
}
else
{
ylodisp <- log(2*seq_len(nbtrials)[-1] - 3)/(log(2))
xlodisp <- c(1, ylodisp - floor(ylodisp))
Umin <- matrix(rep.int(lower, nbtrials), nrow = nbtrials, ncol = nbpar, byrow=T)
etendue <- upper - lower
return(Umin + matrix(rep.int(etendue, nbtrials), nrow = nbtrials, ncol = nbpar, byrow=T)*
sapply(seq_len(nbpar), function (x) sample(xlodisp, nbtrials, replace = FALSE)))
}
}
if (method == "runifbase")
{
return (callrunifmat(nbtrials = nbtrials, nbpar = nbpar, lower = lower, upper = upper))
}
if (method == "runifantithetics")
{
if(floor(nbtrials) != nbtrials)
{
stop("To achieve this, nbtrials should be an even number")
}
return(callrunifmatantith(nbtrials = nbtrials, nbpar = nbpar, lower = lower, upper = upper))
}
if (method == "sobol")
{
if (nbpar == 1)
{
return(lower + (upper - lower) * randtoolbox::sobol(n = nbtrials, dim = nbpar, scrambling = 2))
}
else
{
Umin <- matrix(rep.int(lower, nbtrials), nrow = nbtrials, ncol=nbpar, byrow=T)
etendue <- upper - lower
return(Umin + matrix(rep.int(etendue, nbtrials), nrow = nbtrials, ncol=nbpar, byrow=T)*randtoolbox::sobol(n = nbtrials, dim = nbpar, scrambling = 2))
}
}
if (method == "torus")
{
if (nbpar == 1)
{
return(lower + (upper - lower) * randtoolbox::torus(n = nbtrials, dim = nbpar))
}
else
{
Umin <- matrix(rep.int(lower, nbtrials), nrow = nbtrials, ncol=nbpar, byrow=T)
etendue <- upper - lower
return(Umin + matrix(rep.int(etendue, nbtrials), nrow = nbtrials, ncol=nbpar, byrow=T)*randtoolbox::torus(n = nbtrials, dim = nbpar))
}
}
if (method == "niederreiterlodisp")
{
return(callxlodisp(nbtrials = nbtrials, lower = lower, upper = upper))
}
}
call_localoptim <- function(start0, objectivefn, gradient = NULL, ..., hessian = NULL,
lower, upper, control = list(), themethod, nbtrials = NULL,
verb=FALSE, cl = NULL)
{
nbpar <- length(lower)
if (missing(cl) || is.null(cl))
{
if (missing(nbtrials) || is.null(nbtrials) || nbtrials <= 0)
{
if (themethod == "L-BFGS-B")
{
if (missing(hessian)) hessian <- FALSE
return(optim(par = start0, fn = objectivefn, gr = gradient, ...,
method = themethod, lower = lower, upper = upper,
control = control, hessian = hessian))
}
if (themethod == "Nelder-Mead")
{
if (missing(hessian)) hessian <- FALSE
return(optim(par = start0, fn = objectivefn, gr = gradient, ...,
method = themethod, control = control, hessian = hessian))
}
if (themethod == "nlminb")
{
return(nlminb(start = start0, objective = objectivefn, gradient = gradient,
hessian = hessian, ..., scale = 1, control = control,
lower = lower, upper = upper))
}
}
else
{
ObjectifNSScour <- 1000000
ObjectifNSSi <- 0
res_cour <- 0
nbpar <- length(lower)
if (length(upper) != nbpar || lower > upper) stop("lower and upper are incorrect")
if (verb == TRUE)
{
verbtrace <- list(par = matrix(NA, nrow=nbtrials, ncol=nbpar),
objective = rep.int(NA, nbtrials),
iteration = rep.int(NA, nbtrials))
}
howFar <- txtProgressBar(min=0,max=nbtrials,style=3)
for(i in seq_len(nbtrials))
{
if(nbpar == 1)
{
start_0 <- start0[i]
} else
{
start_0 <- start0[i,]
}
if (themethod == "L-BFGS-B")
{
if (missing(hessian)) hessian <- FALSE
res_cour <- suppressWarnings(try(optim(par = start_0, fn = objectivefn, gr = gradient, ...,
method = themethod, lower = lower, upper = upper,
control = control, hessian = hessian), silent = TRUE))
try(ObjectifNSSi <- res_cour$value, silent = TRUE)
}
if (themethod == "Nelder-Mead")
{
if (missing(hessian)) hessian <- FALSE
res_cour <- suppressWarnings(try(optim(par = start_0, fn = objectivefn, gr = gradient, ...,
method = themethod,
control = control, hessian = hessian), silent = TRUE))
ObjectifNSSi <- res_cour$value
}
if (themethod == "nlminb")
{
res_cour <- suppressWarnings(try(nlminb(start = start_0, objective = objectivefn, gradient = gradient, hessian = hessian,
..., scale = 1, control = control, lower = lower, upper = upper), silent = TRUE))
ObjectifNSSi <- res_cour$objective
}
if(ObjectifNSSi < ObjectifNSScour)
{
res <- res_cour
ObjectifNSScour <- ObjectifNSSi
if (verb == TRUE)
{
verbtrace$par[i, ] <- res_cour$par
verbtrace$objective[i] <- ObjectifNSSi
verbtrace$iteration[i] <- i
}
}
setTxtProgressBar(howFar, value=i)
}
close(howFar)
if (verb == TRUE)
{
booltrace <- !is.na(verbtrace$iteration)
if (nbpar == 1)
{startingparams <- start0[booltrace]}
else startingparams <- start0[booltrace,]
foundparams <- verbtrace$par[booltrace, ]
objective_val <- verbtrace$objective[booltrace]
iteration_no <- verbtrace$iteration[booltrace]
return(list(res=res,
iteration_no = iteration_no,
startingparams_sequence = startingparams,
foundparams_sequence = foundparams,
objective_val_sequence = objective_val))
} else return(res)
}
}
else
{
force(objectivefn)
if (!missing(gradient)) force(gradient)
if (length(upper) != nbpar || lower > upper) stop("lower and upper are incorrect")
cl <- makeCluster(cl, type="SOCK")
.env <- environment()
ldots <- list(...)
clusterExport(cl, "ldots", envir = .env)
if (themethod == "L-BFGS-B")
{
if (missing(hessian)) hessian <- FALSE
if (nbpar <= 1)
{
objectivemin <- suppressWarnings(parSapply(cl, start0, function (x) try(optim(x, fn = objectivefn, gr = gradient, ...,
method = themethod,
lower = lower, upper = upper,
control = control, hessian = hessian), silent=TRUE)))
on.exit(stopCluster(cl))
packageStartupMessage(" Done !")
return(objectivemin[ , which.min(objectivemin[2,])])
}
else
{
objectivemin <- suppressWarnings(parRapply(cl, start0, function (x) try(optim(x, fn = objectivefn, gr = gradient, ...,
method = themethod,
lower = lower, upper = upper,
control = control, hessian = hessian), silent=TRUE)))
on.exit(stopCluster(cl))
packageStartupMessage(" Done !")
return(objectivemin[[unique(which.min(parSapply(cl, objectivemin, '[[', "value")))]])
}
}
if (themethod == "Nelder-Mead")
{
if (missing(hessian)) hessian <- FALSE
if (nbpar <= 1)
{
objectivemin <- suppressWarnings(parSapply(cl, start0, function (x) try(optim(x, fn = objectivefn, gr = gradient, ...,
method = themethod,
control = control, hessian = hessian), silent=TRUE)))
on.exit(stopCluster(cl))
packageStartupMessage(" Done !")
return(objectivemin[ , which.min(objectivemin[2,])])
}
else
{
objectivemin <- suppressWarnings(parRapply(cl, start0, function (x) try(optim(x, fn = objectivefn, gr = gradient, ...,
method = themethod,
control = control, hessian = hessian), silent=TRUE)))
on.exit(stopCluster(cl))
packageStartupMessage(" Done !")
return(objectivemin[[unique(which.min(parSapply(cl, objectivemin, '[[', "value")))]])
}
}
if (themethod == "nlminb")
{
if (!missing(hessian)) force(hessian)
if (nbpar <= 1)
{
objectivemin <- suppressWarnings(parSapply(cl, start0, function (x) try(nlminb(x, objective = objectivefn,
gradient = gradient, hessian = hessian, ...
, scale = 1, control = control,
lower = lower, upper = upper), silent=TRUE)))
on.exit(stopCluster(cl))
packageStartupMessage(" Done !")
return(objectivemin[ , which.min(objectivemin[2,])])
}
else
{
objectivemin <- suppressWarnings(parRapply(cl, start0, function (x) try(nlminb(x, objective = objectivefn,
gradient = gradient, hessian = hessian, ...
, scale = 1, control = control,
lower = lower, upper = upper), silent=TRUE)))
on.exit(stopCluster(cl))
packageStartupMessage(" Done !")
return(objectivemin[[unique(which.min(parSapply(cl, objectivemin, '[[', "objective")))]])
}
}
}
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.