# R/curve.R In curvir: Specify Reserve Demand Curves

#### Documented in curvecurveoptcurvepredinvcurvepredict.curvir

#' Reserve demand curve
#'
#' Fits the reserve demand curve between excess reserves and normalised rates
#'
#' @param y A vector of normalised interest rates.
#' @param x A matrix of explanatory variables. Excess reserve must be the first input.Additional regressor follow (optional).
#' @param type The type of the reserve demand curve. This can be any of \code{logistic}, \code{redLogistic}, \code{fixLogistic}, \code{doubleExp}, \code{exponential}, \code{fixExponential}, \code{arctan}, \code{linear}. See details in \code{\link{curve}}
#' @param dummy Optional input to signify a regime change (vertical shifts in the curve). Must be a vector of equal length to the rows of \code{x}. If not needed use \code{NULL}.
#' @param q Target interval. This is a scalar below 1, for example 0.9 is the 90\% interval. If \code{NULL} then no quantiles are estimated.
#'
#' @note An additional column for the constant is automatically generated, unless requested otherwise.
#'
#' @return Returns a model of class \code{curvir}. This includes
#' \itemize{
#'  \item \code{type} the type of the curve.
#'  \item \code{constant} a logical indicating the use of a constant.
#'  \item \code{w} a list including: \code{mean} the curve parameters for the mean of the curve, \code{upper} and \code{lower} the parameters for the curve at the upper and lower intervals.
#'  \item \code{data} a list including the \code{y}, \code{x}, and \code{dummy} used for the fitting of the curve.
#'  \item \code{mse} the MSE from the fitting of the curve (the mean only).
#'  \item \code{q} the interval used in the fitting of the curve.
#' }
#'
#' @author Nikolaos Kourentzes, \email{nikolaos@kourentzes.com}
#'
#'
#' @references Chen, Z., Kourentzes, N., & Veyrune, R. (2023). \href{https://www.imf.org/en/Publications/WP/Issues/2023/09/01/Modeling-the-Reserve-Demand-to-Facilitate-Central-Bank-Operations-538754}{Modeling the Reserve Demand to Facilitate Central Bank Operations.} IMF Working Papers, 2023(179).
#'
#' @details For a description of the parametric curves, see the provided reference. Below we list their functions:
#' \itemize{
#'  \item \code{logisitc} (Logistic) \deqn{r_i = \alpha + \kappa / (1 - \beta e^{g(\bm{C}_i)}) + \varepsilon_i}
#'  \item \code{redLogistic} (Reduced logistic) \deqn{r_i = \alpha + 1 / (1 - \beta e^{g(\bm{C}_i)}) + \varepsilon_i}
#'  \item \code{fixLogistic} (Fixed logistic) \deqn{r_i = \alpha + 1 / (1 - e^{g(\bm{C}_i)}) + \varepsilon_i}
#'  \item \code{doubleExp} (Double exponential) \deqn{r_i = \alpha + \beta e^{\rho e^{g(\bm{C}_i)}} + \varepsilon_i}
#'  \item \code{exponential} (Exponential) \deqn{r_i = \alpha + \beta e^{g(\bm{C}_i)} + \varepsilon_i}
#'  \item \code{fixExponential} (Fixed exponential) \deqn{r_i = \beta e^{g(\bm{C}_i)} + \varepsilon_i}
#'  \item \code{arctan} (Arctangent) \deqn{r_i = \alpha + \beta \arctan ( g(\bm{C}_i))  + \varepsilon_i}
#'  \item \code{linear} (Linear) \deqn{r_i = g(\bm{C}_i) + \varepsilon_i}
#' }
#' And \eqn{g(\bm{C}) = c + \bm{C} w_g}, where \eqn{\alpha}, \eqn{\beta}, \eqn{\kappa}, \eqn{\rho} are curve parameters,
#' \eqn{c} is a constant togglable by \code{constant}, \eqn{\bm{C}} are the regressors including the excess reserves. \eqn{w_g} their coefficients, and finally \eqn{\varepsilon_i} is the error term of the curve.
#'
#' @examples
#' \dontshow{
#'   rate <- head(ecb$rate,10) #' x <- ecb$x[1:10,1,drop=FALSE]
#'   curve(x,rate,rep=1,type="fixExponential")
#' }
#' \donttest{
#'   # Use ECB example data
#'   rate <- ecb$rate #' x <- ecb$x[,1,drop=FALSE]
#'   curve(x,rate)
#'
#'   # An arctangent curve
#'   curve(x,rate,type="arctan")
#'  }
#'
#' @export curve

curve <- function(x,y,type="logistic",dummy=NULL,q=NULL,...){

type <- match.arg(type,c("logistic","redLogistic","fixLogistic",
"doubleExp","exponential","fixExponential",
"arctan",
"linear"))

# Optimise curve
wOpt <- curveopt(x=x,y=y,type=type,dummy=dummy,q=NULL,...)

# Extract arguments from ellipsis and if missing use defaults from curveopt()
elarguments <- list(...)
if ((length(elarguments)>0) && any(grepl("constant",names(elarguments)))){
constant <- elarguments$constant } else { constant <- eval(formals(curveopt)$constant)[1]
}

# Estimate quantiles
if (!is.null(q)){
# Get fitted values
if (constant){
xhat <- cbind(1,x)
} else {
xhat <- x
}
yhat <- curvepred(x=xhat,w=wOpt$w,type=type,dummy=dummy) # Optimise # Restricted to not cross yhat p <- (1-q)/2 wU <- curveopt(x,y,type=type,q=1-p,winit=wOpt$w,yhat=yhat,dummy=dummy,constant=constant,...)$w wL <- curveopt(x,y,type=type,q=p,winit=wOpt$w,yhat=yhat,dummy=dummy,constant=constant,...)$w } else { wU <- wL <- NULL } # Save curve model data = list(y=y,x=x,dummy=dummy) w <- list(mean=wOpt$w,upper=wU,lower=wL)
model <- list(type=type,constant=constant,w=w,data=data,mse=wOpt$mse,q=q) model <- structure(model,class="curvir") return(model) } #' Predict method for curvir reserve demand models #' #' Predicted values based on curvir model object #' #' @param object A model fit with \code{\link{curve}}. #' @param newdata New input data organised as the x matrix in \code{\link{curve}}. If \code{NULL} then the data used to fit the model is re-used. #' @param newdummy New input dummy organised as the dummy vector in \code{\link{curve}}. If \code{NULL} then the dummy used to fit the model is re-used. #' @param ... Further arguments (unused) #' #' @return Returns a matrix of predicted values. If the model has estimates for intervals then it will provide upper and lower intervals. #' #' @inherit curve references author #' @seealso \code{\link{curve}}. #' #' @examples #' \dontshow{ #' rate <- head(ecb$rate,10)
#'   x <- ecb$x[1:10,1,drop=FALSE] #' fit <- curve(x,rate,rep=1,type="fixExponential") #' predict(fit) #' } #' \donttest{ #' # Use ECB example data #' rate <- ecb$rate
#'   x <- ecb$x[,1,drop=FALSE] #' fit <- curve(x,rate) #' predict(fit) #' #' # An example with new data #' predict(fit,newdata=tail(x)) #' } #' #' @describeIn predict Predicted values for parametric curves #' @export #' @method predict curvir predict.curvir <- function(object, newdata=NULL, newdummy=NULL,...){ # Use old data if nothing new is provided if (is.null(newdata)){ newdata <- object$data$x } if (is.null(newdummy)){ newdummy <- object$data$dummy } # Add constant as needed if (object$constant){
newdata <- cbind(1,newdata)
}

if (!is.null(newdummy)){
if (length(newdummy) != nrow(newdata)){
stop("Length of newdummy must match the rows of newdata.")
}
}

yhatMean <- curvepred(newdata,object$w$mean,object$type,dummy=newdummy) colnames(yhatMean) <- "mean" # Interval if (!is.null(object$w$upper)){ yhatUpper <- curvepred(newdata,object$w$upper,object$type,dummy=newdummy)
colnames(yhatUpper) <- paste0("upper ",object$q*100,"%") } else { yhatUpper <- NULL } if (!is.null(object$w$lowe)){ yhatLower <- curvepred(newdata,object$w$lower,object$type,dummy=newdummy)
colnames(yhatLower) <- paste0("lower ",object$q*100,"%") } else { yhatLower <- NULL } return(cbind(yhatMean,yhatLower,yhatUpper)) } #' Reserve demand curve predicted values #' #' Provides the predicted values for the reserve demand curve of choice. For general use prefer the #' predict() function, which handles the constant internally. #' #' @param x A matrix with the inputs. If there is a constant in the estimated curve, then the first column in \code{x} should be the constant. Next column should be the excess reserves. Additional regressor follow (optional). #' @param w The vector of weights for the desired curve. Estimated using the \code{\link{curveopt}} function. #' @inheritParams curve #' #' @return Returns a vector of the predicted values. #' #' @seealso \code{\link{curve}}, and \code{\link{curveopt}}. #' @inherit curve references details author #' @export curvepred # Note: Constant must be the first column in x, and excess reserve the next curvepred <- function(x,w,type="logistic",dummy=NULL){ type <- match.arg(type,c("logistic","redLogistic","fixLogistic", "doubleExp","exponential","fixExponential", "arctan", "linear")) p <- ncol(x) b <- tail(w,p) g <- gX(x,b) yhat <- switch(type, logistic = logistic(w,g,dummy), redLogistic = redLogistic(w,g,dummy), fixLogistic = fixLogistic(w,g,dummy), doubleExp = doubleExp(w,g,dummy), exponential = exponential(w,g,dummy), fixExponential = fixExponential(w,g,dummy), arctan = arctan(w,g,dummy), linear = linear(w,g,dummy) ) return(yhat) } #' Optimise curve parameters #' #' Finds optimal curve parameters. #' #' @inheritParams curve #' @param constant A logical (\code{TRUE} or \code{FALSE}) whether to include a constant or not. #' @param reps Number of repetitions for the particle swarm optimisation. #' @param sign A vector of equal length to the number of additional regressors in \code{x} (excluding the constant (if used) and the excess reserves) of positive and negative values (any) that will be used to obtain signs to restrict the values of the estimated parameters. Use \code{NULL} for no restrictions. #' @param q The desired quantile to optimise for. Use \code{NULL} to get the conditional expectation. #' @param winit A vector of initial values for the optimisation. This will also carry over to sign restrictions if \code{sameSign==TRUE}. #' @param yhat Useful when estimating quantiles. Supply here the predicted values for the conditional expectation to add restrictions for the quantiles to not cross the conditional expectation. Use \code{NULL} to not add any restrictions. #' @param wsel Use the minimum error set of parameters (\code{select}) or the combination of the pool parameters using the heuristic in Kourentzes et al., (2019) (\code{combine}). #' @param sameSign Used if \code{winit != NULL} to take any sign restrictions from \code{winit}. #' #' @return Returns a list of #' \itemize{ #' \item \code{w} The optimal parameters #' \item \code{mse} The Mean Squared Error of the fitted curve. #' } #' #' @seealso \code{\link{curve}}, and \code{\link{curvepred}}. #' @inherit curve author #' @references \itemize{ #' \item Chen, Z., Kourentzes, N., & Veyrune, R. (2023). \href{https://www.imf.org/en/Publications/WP/Issues/2023/09/01/Modeling-the-Reserve-Demand-to-Facilitate-Central-Bank-Operations-538754}{Modeling the Reserve Demand to Facilitate Central Bank Operations.} IMF Working Papers, 2023(179). #' \item Kourentzes, N., Barrow, D., & Petropoulos, F. (2019). Another look at forecast selection and combination: Evidence from forecast pooling. International Journal of Production Economics, 209, 226-235. #' } #' #' @examples #' \dontshow{ #' rate <- head(ecb$rate,10)
#'   x <- ecb$x[1:10,1,drop=FALSE] #' curveopt(x,rate,reps=1,type="fixExponential") #' } #' \donttest{ #' # Use ECB example data #' rate <- ecb$rate
#'   x <- ecb$x[,1,drop=FALSE] #' curveopt(x,rate) #' } #' #' @export curveopt # Function to optimise curve parameters curveopt <- function(x,y,type="logistic",constant=c(TRUE,FALSE), reps=3,sign=NULL,q=NULL,winit=NULL,yhat=NULL, wsel=c("select","combine"),dummy=NULL, sameSign=c(TRUE,FALSE)){ constant <- constant[1] sameSign <- sameSign[1] wsel <- match.arg(wsel,c("select","combine")) # Optimise coefficients # Clean data idx <- !(apply(x,1,function(x){any(is.na(x))}) | is.na(y)) xClean <- x[idx,,drop=FALSE] xClean <- as.matrix(cbind(xClean)) yClean <- y[idx] if (!is.null(dummy)){ dummy <- dummy[idx] } # # Order data # idx <- order(xClean[,1]) # xClean <- xClean[idx,,drop=FALSE] # yClean <- yClean[idx] # Normalise variance - needed for better optimisation sdx <- apply(xClean,2,sd) xCleanSc <- xClean / matrix(rep(sdx,nrow(xClean)),ncol=ncol(xClean),byrow=TRUE) # xClean[,1] <- xClean[,1] / sdx[1] # Introduce constant if (constant){ X <- cbind(1,xCleanSc) } p <- ncol(X) # p includes the constant eta <- as.numeric(!is.null(dummy)) if (is.null(winit)){ winit <- switch(type, logistic = c(0,1,0,rep(-0.1,p+eta)), redLogistic = c(0,1,rep(-0.1,p+eta)), fixLogistic = c(0,rep(-0.1,p+eta)), doubleExp = c(0.2,0,-0.1,1,rep(0,p-1+eta)), exponential = c(0.2,0,1,rep(0,p-1+eta)), fixExponential = c(0.2,1,rep(0,p-1+eta)), arctan = c(1,-1,rep(0,p+eta)), linear = rep(0,p+eta) ) } else { winit <- winit # If a winit is provided and sameSign is TRUE than retain # signs for the X's if (sameSign){ sign <- tail(sign(winit),p-constant) } } # Restrict for signs of explanatory variables upper <- rep(20,length(winit)) lower <- rep(-20,length(winit)) if (!is.null(sign)){ upper[length(winit) - (ncol(x) + constant + eta) + 1 + which(sign < 0)] <- 0 lower[length(winit) - (ncol(x) + constant + eta) + 1 + which(sign > 0)] <- 0 } # Pre-optmise using Nelder-Meade opt1 <- optim(winit,loss,y=yClean,x=X,type=type,q=q,yhat=yhat,dummy=dummy) # Override the solution of pre-optimisation if it violates the bounds opt1$par[opt1$par>upper] <- upper[opt1$par>upper]
opt1$par[opt1$par<lower] <- lower[opt1$par<lower] # We will optimise multiple times to reduce the stochasticity of the results wOptAll <- array(NA,c(reps,length(winit))) fnAll <- vector("numeric",reps) for (r in 1:reps){ # Optimise using particle swarm psopts <- list(maxit=4000,type="SPSO2007",fnscale=1,s=20) opts <- psoptim(opt1$par, loss, upper=upper, lower=lower, control=psopts,
y=yClean ,x=X ,type=type, q=q, yhat=yhat, dummy=dummy)

wOptAll[r,] <- opts$par fnAll[r] <- opts$value
}

repsRnk <- order(fnAll)
wOptAll <- wOptAll[repsRnk,,drop=FALSE]
fnAll <- fnAll[repsRnk]

if (wsel == "combine"){
# Use islands to select what to combine
if (reps != 1){

# Get correlation from min MSE solution
crit <- -1*(cor(t(wOptAll))[1,])
# crit <- c(0,diff(fnAll))

# Re-order based on critical value
critRnk <- order(crit)
crit <- crit[critRnk]
wOptAll <- wOptAll[critRnk,,drop=FALSE]
fnAll <- fnAll[critRnk]

# Combination islands
critD <- c(0,diff(crit))
# Create a list of increasing elements
crt <- lapply(1:length(critD),function(x){critD[1:x]})
s <- 1.5 # Factor to scale IQR
th <- sapply(crt,function(x){s*IQR(x)+quantile(x,0.75)})
th[1] <- NA
# Find where the threshold is crossed
cKeep <- which(critD>th)
if (length(cKeep)==0){
cKeep <- reps
} else {
cKeep <- cKeep[1]-1
}
cKeep <- max(cKeep,1)

# Override for low correlation
cKeep <- min(cKeep,tail(which(cor(t(wOptAll))[1,]>=0.8),1))

# plot(X[,1+constant],yClean)
# cmpT <- colorRampPalette(RColorBrewer::brewer.pal(9,"Reds")[3:8])(reps)
# for (i in 1:reps){
#   ord <- order(X[,1+constant])
#   temp <- curvepred(X,wOptAll[i,],type=type)
#   lines(X[ord,1+constant],temp[ord],col=cmpT[i])
# }
# legend("topright",paste0("rep",1:reps),col=cmpT,bty="n",cex=0.7,lty=1)

# plot(log(critD),type="o",pch=20)
# lines(log(th),col="red")
# abline(v=cKeep)

# Do a weighted combination on fnAll (MSE)
keepIdx <- c(rep(TRUE,cKeep),rep(FALSE,reps-cKeep))
temp <- wOptAll[keepIdx,,drop=FALSE]
# wOpt <- as.numeric(rbind((1/(crit[keepIdx]+1.1))/sum(1/(crit[keepIdx]+1.1))) %*% temp)
wOpt <- as.numeric(rbind(1/fnAll[keepIdx]/sum((1/fnAll[keepIdx]))) %*% temp)
} else {
# cKeep <- 1
wOpt <- wOptAll
}
} else {
# Select only top performing set of coefficients
wOpt <- wOptAll[which.min(fnAll),]
}

# Re-scale parameters
wOpt[(length(wOpt)-ncol(xClean)+1):length(wOpt)] <- wOpt[(length(wOpt)-ncol(xClean)+1):length(wOpt)] / sdx
# wOpt[(length(wOpt)-ncol(xClean)+1)] <- wOpt[(length(wOpt)-ncol(xClean)+1)] / sdx[1]

if (constant){
xhat <- cbind(1,xClean)
}

# ord <- order(xhat[,1+constant])
# plot(xhat[ord,1+constant],yClean[ord])
# lines(xhat[ord,1+constant],
#       curvepred(xhat[ord,,drop=FALSE],wOpt,type),col="red")

# Calculate in-sample MSE on original scale
mse <- mean((yClean - curvepred(xhat[,,drop=FALSE],wOpt,type,dummy))^2)

# Name parameters
wOpt <- namePar(wOpt,type,xClean,constant,dummy)

return(list(w=wOpt,mse=mse))

}

#' Calculate the inverse curve prediction
#'
#' Calculate the predicted reserves given some rate, i.e., calculate the prediction of the inverse curve.
#'
#' @param object A model fit with \code{\link{curve}}.
#' @param ynew The input rate. If \code{NULL} this corresponds to the values from \code{predict(object)}.
#' @param xnew The values for the additional regressors that were used in the curve fit. Must be a matrix, ordered (columns) as they were input in the fitting of the curve. The constant is dealt with automatically. Do not input the excess reserves. If \code{NULL} this is picked up from the data used to fit the curve.
#' @param dummynew The values for the indicator, if one was used in the fitting of the curve. If \code{NULL} then the data used in the fitting of the model are used.
#' @param warn A logical (\code{TRUE} or \code{FALSE}) to issue a warning if the resulting values are more than 10\% away from the min-max of the excess reserves used to estimate the curve.
#'
#' @return Returns a vector of values of the predicted reserves
#'
#' @inherit curve author references
#'
#' @examples
#' \dontshow{
#'   rate <- ecb$rate #' x <- ecb$x[,1,drop=FALSE]
#'   fit <- curve(x,rate,rep=1,type="fixExponential")
#'   invcurve(fit)
#' }
#' \donttest{
#'   # Use ECB example data
#'   rate <- ecb$rate #' x <- ecb$x[,1,drop=FALSE]
#'   fit <- curve(x,rate,type="logistic")
#'   invcurve(fit)
#'
#'   # Use a different input rate
#'   invcurve(fit,ynew=0.1)
#' }
#'
#' @export invcurve

invcurve <- function(object,ynew=NULL,xnew=NULL,dummynew=NULL,warn=c(TRUE,FALSE)){

type <- object$type w <- object$w$mean constant <- object$constant
warn <- warn[1]

# Keep the mean prediction
yhat <- predict(object)[,1,drop=FALSE]

if (is.null(ynew)){
# Use data in the curve object
ynew <- yhat
if (is.null(xnew)){
xnew <- object$data$x
# Remove the excess reserves
xnew <- xnew[,-1,drop=FALSE]
}
if (is.null(ynew)){
dummynew <- object$data$dummy
}
} else {
# Data in ynew is external, make sure that the
# correct inputs are provided
p <- ncol(object$data$x)-1 # remove the excess reserves
if (!is.null(xnew)){
if (ncol(xnew) != p){
stop("xnew must be a matrix with the same number of explanatory variables as the matrix used to fit the curve. Do not input the excess reserve here.")
}
if (length(ynew) != nrow(xnew)){
stop("xnew must have equal number of rows as the values in ynew.")
}
}
# Check the new dummy
if (!is.null(object$data$dummy)){
if (is.null(dummynew)){
stop("The curve has been fitted with a dummy. Please input a dummy information.")
}
if (length(ynew) != length(dummynew)){
stop("dummy must have the same length as ynew.")
}
}
}

invy <- switch(type,
logistic = invlogistic(w,ynew),
redLogistic = invredLogistic(w,ynew),
fixLogistic = invfixLogistic(w,ynew),
doubleExp = invdoubleExp(w,ynew),
exponential = invexponential(w,ynew),
fixExponential = invfixExponential(w,ynew),
arctan = invarctan(w,ynew),
linear = invlinear(w,ynew)
)

# Inverse the g(C) part for given explanatory values
# First construct the value of all X's without the reserves
# Obtain coefficients
if (!is.null(xnew)){
p <- ncol(xnew) + constant
} else {
p <- as.numeric(constant)
}
b <- tail(w,p+1) # include +1 for the excess reserves
if (constant){
xnew <- cbind(1,xnew)
}
# xalpha is the "constant" that goes in the inverse of gX
# This incorporates the effect of all explanatories, the constant, and the dummy
# xbeta is the coefficient for the excess reserves
xalpha <- xnew %*% b[-(1+constant)]
if (!is.null(object$data$dummy)){
eta <- w[4]
xalpha = eta * dummynew
}
xbeta <- b[(1+constant)]
invyy <- invgX(invy,xalpha,xbeta)

if (any(invyy < min(object$data$x[,1]) - 0.1*diff(range(object$data$x[,1]))) |
any(invyy > max(object$data$x[,1]) + 0.1*diff(range(object$data$x[,1]))) &&
warn){
warning("Predicted value(s) is more than 10% away from the min/max of the data used for fitting the curve.")
}

return(invyy)

}

## Internal functions

#' @export
#' @method print curvir
print.curvir <- function(x, ...){
type <- x$type writeLines(paste0(toupper(substr(x$type,1,1)), tolower(substr(x$type,2,nchar(x$type))), " type reserve demand curve"))
if (ncol(x$data$x)-1>0){
writeLines(paste0(ncol(x$data$x)-1, " additional regressors included."))
}
if (!is.null(x$dummy)){ writeLines("Shift indicator included.") } writeLines("") writeLines("Curve parameters:") print(x$w$mean) if (!is.null(x$q)){
writeLines("")
writeLines(paste0("Model has estimates for ",round(x\$q*100,0),"% intervals."))
}

}

# Name curve parameters
namePar <- function(wOpt,type="logistic",x=NULL,constant=c(TRUE,FALSE),dummy=NULL){

constant <- constant[1]
type <- match.arg(type,c("logistic","redLogistic","fixLogistic",
"doubleExp","exponential","fixExponential",
"arctan",
"linear"))

wNm <- switch(type,
logistic = c("alpha","beta","kappa"),
redLogistic = c("alpha","beta"),
fixLogistic = c("alpha"),
doubleExp = c("alpha","beta","rho"),
exponential = c("alpha","beta"),
fixExponential = c("beta"),
arctan = c("alpha","beta"),
linear = NULL
)
# If there is a dummy name eta
if (!is.null(dummy)){
wNm <- c(wNm,"eta")
}
# If there is a constant, name it
if (constant){
wNm <- c(wNm,"constant")
}
# Names of X's
if (!is.null(x)){
wNm <- c(wNm,colnames(x))
} else {
xNm <- paste0("X",1:(length(wOpt)-length(wNm)))
xNm[1] <- "Excess"
wNm <- c(wNm,xNm)
}

names(wOpt) <- wNm
return(wOpt)

}

# Loss function for optimising curve parameters
loss <- function(w,y,x,type,q=NULL,yhat=NULL,dummy=NULL){

# w contains the parameters
# y is the target
# x is the explanatories as a matrix
# q is used to optimise for a desired quantile

yfit <- curvepred(x,w,type,dummy)
e <- y - yfit

if (is.null(q)){
e <- sum(e*e)
} else {
b <- rep(q,length(y))
b[y < yfit] <- 1 - q
e <- sum(b * abs(y - yfit))
}

# Make quantiles to not cross expectation
if (!is.null(q) && !is.null(yhat)){
d <- yfit - yhat
if (any(sign(q-0.5) * d < 0)){
e <- (e+1)^5
}
}

# e <- sum(abs(e)) + diff(range(e))

# # Ensure monotonicity
# if (!any(is.na(yfit)) & !any(is.infinite(yfit))){
#   if (any(sign(diff(yfit))<=0)){
#     e <- (e+1)^5
#   }
# }

# if (is.na(e)){
#   e <- 10^10
# }

return(e)

}

# Various curves follow, call them through curvepred()
logistic <- function(w,g,dummy){

alpha <- w[1]
beta <- w[2]
kappa <- w[3]

yhat <- alpha + kappa/(1-beta*exp(g))

if (!is.null(dummy)){
eta <- w[4]
yhat <- yhat + eta*dummy
}

return(yhat)

}

# Estimate g(X) - the exponent part
gX <- function(x,b){
# If there is a constant, add it to X before it is fed here
g <- x %*% b
return(g)
}

# Various curves follow, call them through curvepred()
logistic <- function(w,g,dummy){

alpha <- w[1]
beta <- w[2]
kappa <- w[3]

yhat <- alpha + kappa/(1-beta*exp(g))

if (!is.null(dummy)){
eta <- w[4]
yhat <- yhat + eta*dummy
}

return(yhat)

}

redLogistic <- function(w,g,dummy){

alpha <- w[1]
beta <- w[2]

yhat <- alpha + 1/(1-beta*exp(g))

if (!is.null(dummy)){
eta <- w[3]
yhat <- yhat + eta*dummy
}

return(yhat)

}

fixLogistic <- function(w,g,dummy){

alpha <- w[1]

yhat <- alpha + 1/(1-exp(g))

if (!is.null(dummy)){
eta <- w[2]
yhat <- yhat + eta*dummy
}

return(yhat)

}

doubleExp <- function(w,g,dummy){

alpha <- w[1]
beta <- w[2]
rho <- w[3]

yhat <- alpha + beta * exp(rho*exp(g))

if (!is.null(dummy)){
eta <- w[4]
yhat <- yhat + eta*dummy
}

return(yhat)

}

exponential <- function(w,g,dummy){

alpha <- w[1]
beta <- w[2]

yhat <- alpha + beta * exp(g)

if (!is.null(dummy)){
eta <- w[3]
yhat <- yhat + eta*dummy
}

return(yhat)

}

fixExponential <- function(w,g,dummy){

beta <- w[1]

yhat <- beta*exp(g)

if (!is.null(dummy)){
eta <- w[2]
yhat <- yhat + eta*dummy
}

return(yhat)

}

linear <- function(w,g,dummy){

yhat <- g

if (!is.null(dummy)){
eta <- w[1]
yhat <- yhat + eta*dummy
}

return(yhat)

}

arctan <- function(w,g,dummy){

alpha <- w[1]
beta <- w[2]

yhat <- alpha + beta*atan(g)

if (!is.null(dummy)){
eta <- w[3]
yhat <- yhat + eta*dummy
}

return(yhat)

}

invlogistic <- function(w,ynew){

alpha <- w[1]
beta <- w[2]
kappa <- w[3]

# Bound by alpha
ynew[ynew<alpha] <- alpha + 1e-18
# Deal with the asymptote
dy <- (alpha - ynew)
dy[dy == 0] <- -1e-18

# Find inverse
yinv <- log(-(-kappa/(dy) - 1)/beta)

return(yinv)

}

# Estimate g(X) - the exponent part
invgX <- function(x,alpha,beta){
# If there is a constant, add it to X before it is fed here
g <- (x - alpha)/beta
return(g)
}

invredLogistic <- function(w,ynew){

alpha <- w[1]
beta <- w[2]

# Bound by alpha
ynew[ynew<alpha] <- alpha + 1e-18
# Deal with the asymptote
dy <- (alpha - ynew)
dy[dy == 0] <- -1e-18
# Find inverse
yinv <- log(-(-1/(dy) - 1)/beta)

return(yinv)

}

invfixLogistic <- function(w,ynew){

alpha <- w[1]

# Bound by alpha
ynew[ynew<alpha] <- alpha + 1e-18
# Deal with the asymptote
dy <- (alpha - ynew)
dy[dy == 0] <- -1e-18
# Find inverse
yinv <- log(1/(dy) + 1)

return(yinv)

}

invdoubleExp <- function(w,ynew){

alpha <- w[1]
beta <- w[2]
rho <- w[3]

# Deal with the asymptote
# Deal with negatives
dy <- (ynew-alpha)/beta
dy[dy <= 0] <- 1e-18
# Deal with negatives again
dy <- log(dy)/rho
dy[dy <= 0] <- 1e-18
# Find inverse
yinv <- log(dy)

return(yinv)

}

invexponential <- function(w,ynew){

alpha <- w[1]
beta <- w[2]

# Deal with the asymptote
dy <- (ynew-alpha)/beta
dy[dy <= 0] <- 1e-18
# Find inverse
yinv <- log(dy)

return(yinv)

}

invfixExponential <- function(w,ynew){

beta <- w[1]
# Find inverse
yinv <- log(ynew/beta)

return(yinv)

}

invlinear <- function(w,ynew){

yinv <- ynew

return(yinv)

}

invarctan <- function(w,ynew){

alpha <- w[1]
beta <- w[2]

yinv <- tan((ynew - alpha)/beta)
return(yinv)

}

# Outputs the number of parameters in a curve
modelSize <- function(x,type,constant=c(TRUE,FALSE),dummy=NULL){

constant <- constant[1]
constant <- as.numeric(constant)
eta <- as.numeric(!is.null(dummy))

p <- switch(type,
logistic = ncol(x)+3+constant+eta,
redLogistic = ncol(x)+2+constant+eta,
fixLogistic = ncol(x)+1+constant+eta,
doubleExp = ncol(x)+3+constant+eta,
exponential = ncol(x)+2+constant+eta,
fixExponential = ncol(x)+1+constant+eta,
arctan = ncol(x)+2+constant+eta,
linear = ncol(x)+constant+eta
)

return(p)

}


## Try the curvir package in your browser

Any scripts or data that you put into this service are public.

curvir documentation built on Nov. 24, 2023, 5:09 p.m.