This is my code i don't understand the error message:
library(rgenoud)
rm(list=ls())
set.seed(666)
#########################################################
# As a first step, it is assumed that all input parameters are independent of
ageing :
#########################################################
InputDim <-20
# Max number of ageings in the inputs
CPIRate <- rep(0.02 , InputDim )
# CPI inflation rate
LossInflaRate <- rep(0.02 , InputDim)
# inflation rate of losses
RiskFreeRate <- rep(0.04, InputDim )
# vector of risk free rate
RiskRate <- rep(0.06, InputDim )
# vector of risk free rate
################ lapse param begin :
LapseType <- rep("linear" , InputDim )
# type of the function used to calculate lapse rate= "power", "logistic",
"linear"
CentralLapseRate <- rep(0.13 , InputDim )
# lapse rate for price factor = 1
ElasticityLeft <- rep( 10, InputDim )
# left elasticity for price factor = 1
ElasticityRight <- rep(10, InputDim )
# right elasticity for price factor = 1
WindowLeft <- rep( 1 , InputDim )
# friction window price factor left bound
WindowRight <- rep( 1, InputDim )
# friction window price factor right bound
MinLapseRate <- rep( 0 , InputDim )
# min lapse rate
MaxLapseRate <- rep( 1 , InputDim )
# min lapse rate
Alpha <- 0.5
#acquisition operator
############### lapse param end
Price <- rep(1 , InputDim )
# Axa price for central market position and price factor = 1
IndivMeanLossAmount <- rep(0.97 , InputDim )
# mean individual loss amount
Commission <- rep(0 , InputDim )
# commision rate
IndividualExp <- rep(0 , InputDim )
# expense amount per risk
#######################################
# Market Expert Model parameter begin :
#######################################
LRUpBound <- 1.01349495023810
# Upper bound of the loss ratio
LRLowBound <- 0.743806254238095
# Lower bound of the loss ratio
TransitCoefMean <- 1.03001118178481
# Mean of the coefficient of transition
TransitCoefStdDev <- 0.0187801387384363
# Standard deviation of the coefficient of transition
#######################################
# Optimization algorithm & Monte-Carlo parameters :
#######################################
LTMeanComputYearNb <- 1000
# Number of year used to simulate LR
NbOfSimScenar <- 1000
# Number Monte-Carlo scenarios
GenAlgoIterNb <- 100
# number of iteration for genetic algorithm
#######################################
# Problem size parameters :
#######################################
FirstYear <-2009
# first year of renewal
LastYear <- 2018
# last year of projected renewal
######################################
# Function giving the minimum x (price factor = 1+x ) :
######################################
Xmin <- function(LapseType, ElasticityLeft, WindowLeft, MinLapseRate
,CentralLapseRate)
{
InputDim <- length(LapseType)
x=numeric(InputDim)
for (i in 1:InputDim)
{
if (LapseType[i]=="power" )
x[i] <- WindowLeft [i]*(MinLapseRate[i]
/CentralLapseRate[i])^(1/ElasticityLeft[i])-1
else if (LapseType[i]=="linear")
x [i]<-
(WindowLeft[i]/ElasticityLeft[i])*(ElasticityLeft[i]-1+MinLapseRate[i]
/CentralLapseRate[i])-1
}
x
}
######################################
# Function giving the maximum x (price factor = 1+x ) :
######################################
Xmax <- function(LapseType, ElasticityLeft, WindowLeft, MaxLapseRate
,CentralLapseRate)
{
InputDim <- length(LapseType)
x=numeric(InputDim)
for (i in 1:InputDim)
{
if (LapseType[i]=="power" )
x[i] <- WindowLeft[i]
*(MaxLapseRate[i]/CentralLapseRate[i])^(1/ElasticityLeft[i])-1
else if (LapseType[i]=="linear")
x[i] <-
(WindowLeft[i]/ElasticityLeft[i])*(ElasticityLeft[i]-1+MaxLapseRate[i]/CentralLapseRate[i])-1
}
x
}
##########################################################
# Function simulating a scenario of projected market loss ratio :
##########################################################
LRScenar <- function(PriorLR,Phase,LRUpBound
,LRLowBound,TransitCoefMean,TransitCoefStdDev,LTMeanComputYearNb)
{
initial <- PriorLR
LRScenar <- NULL
phase2 <- Phase
SigmaLog <- (log(1+(TransitCoefStdDev /TransitCoefMean)^2))^0.5
MuLog <- log(TransitCoefMean)-SigmaLog^2/2
for (i in 1:LTMeanComputYearNb)
{
kk <- rlnorm(1, meanlog = MuLog, sdlog = SigmaLog)
if (phase2!=1) kk <- (1/kk)
x_1 <- initial*kk
initial <- x_1
if (((x_1 > LRUpBound )&(phase2==1))|((x_1<LRLowBound)&(phase2==-1))) phase2
<- (-phase2)
LRScenar <- c(LRScenar,x_1)
}
LRScenar
}
################################################
# Function giving the lapse rate :
################################################
LapseRate <- function(x, LapseType,CentralLapseRate,
ElasticityLeft, ElasticityRight,WindowLeft, WindowRight,
MinLapseRate,MaxLapseRate )
{
PriceFactor <- 1+x
if (LapseType=="power")
{
if(PriceFactor<WindowLeft)
r <- CentralLapseRate*(PriceFactor/WindowLeft)^ElasticityLeft
else if((PriceFactor<=WindowRight)&(PriceFactor>=WindowLeft ))
r <- CentralLapseRate
else
r <-CentralLapseRate*(PriceFactor/WindowRight)^ElasticityRight
}
else if (LapseType=="linear")
{
if(PriceFactor<WindowLeft)
r <- CentralLapseRate*(1+ElasticityLeft/WindowLeft*(PriceFactor-WindowLeft))
else if((PriceFactor<=WindowRight)&(PriceFactor>=WindowLeft ))
r <- CentralLapseRate
else
r <-
CentralLapseRate*(1+ElasticityRight/WindowRight*(PriceFactor-WindowRight))
}
min( max( r, MinLapseRate ) , MaxLapseRate)
}
######################################################
# Function calculating a matrix of market price (excluding inflation)
scenarios :
######################################################
MarketPriceMatrix <- function(NbOfFutureSimYear ,NbOfSimScenar,Phase,LRUpBound
,LRLowBound,TransitCoefMean,
TransitCoefStdDev,PriorMarketPricePosition,MarketPricePhase,MarketPriceLTMean)
{
LRPhase <- -MarketPricePhase
MarketPriceMatrix <- matrix(numeric(NbOfFutureSimYear*NbOfSimScenar),
nrow=NbOfSimScenar,ncol=NbOfFutureSimYear)
PriorLR <- 1 / (MarketPriceLTMean * PriorMarketPricePosition)
for (i in 1:NbOfSimScenar)
MarketPriceMatrix[i,] <- 1/ LRScenar (PriorLR,LRPhase,LRUpBound
,LRLowBound,TransitCoefMean,
TransitCoefStdDev,NbOfFutureSimYear)/MarketPriceLTMean
MarketPriceMatrix
}
#############################################
# Function calculating a scenario of results :
#############################################
ResultScenar <-
function(FirstYear,OptimYear,Price,CumulRetentionCoef,Commission,IndividualExp,NbOfFutureSimYear,
IndivMeanLossAmount,CPIRate,LossInflaRate,Alpha)
{
CumulLossInflaRate <- cumprod(1+LossInflaRate)
CumulCPIRate <- cumprod(1+CPIRate)
y <- numeric(NbOfFutureSimYear)
delta <- OptimYear - FirstYear
if (OptimYear==FirstYear)
{
y[1] <CumulRetentionCoef[1]*( (Price[1]*(1-Commission[delta
+1])-IndivMeanLossAmount[delta +1])*CumulLossInflaRate[delta +1]
- IndividualExp[delta +1]*CumulCPIRate[delta +1])
}
for (i in 2: NbOfFutureSimYear)
{
k <- delta + i
y[i] <- (1-Alpha)*CumulRetentionCoef[i-1]*(
(Price[i-1]*(1-Commission[k-1])-IndivMeanLossAmount[k-1])*CumulLossInflaRate[k-1]
- IndividualExp[k]*CumulCPIRate[k-1] )
+ Alpha*CumulRetentionCoef[i]*(
(Price[i]*(1-Commission[k])-IndivMeanLossAmount[k])*CumulLossInflaRate[k]
- IndividualExp[k]*CumulCPIRate[k] )
}
if (OptimYear!=FirstYear)
{
for (i in 1: (NbOfFutureSimYear+1))
{
k <- delta + i
y[i] <- (1-Alpha)*CumulRetentionCoef[i]*(
(Price[i]*(1-Commission[k-1])-IndivMeanLossAmount[k-1])*CumulLossInflaRate[k-1]
- IndividualExp[k]*CumulCPIRate[k-1] )
+ Alpha*CumulRetentionCoef[i+1]*(
(Price[i+1]*(1-Commission[k])-IndivMeanLossAmount[k])*CumulLossInflaRate[k]
- IndividualExp[k]*CumulCPIRate[k] )
}
}
y
}
##############################################
### Function calculating a scenario of results :
###############################################
ResultActuSum <- function(FirstYear,OptimYear,x, NbOfFutureSimYear,
CentralLapseRate, ElasticityLeft,WindowLeft, WindowRight,
ElasticityRight,
MarketPriceMean,Price,Commission,IndividualExp,IndivMeanLossAmount,
CPIRate,LossInflaRate,Beta,Alpha)
{
R <- numeric(NbOfFutureSimYear+1)
Price2 <- numeric(NbOfFutureSimYear+1)
CumulRetentionCoef <- numeric(NbOfFutureSimYear+1)
Result <- numeric(NbOfFutureSimYear+1)
delta <- OptimYear - FirstYear
for (i in 1:(NbOfFutureSimYear+1))
{
k <- delta + i
R[i] <- LapseRate(x[i], LapseType[k],CentralLapseRate[k],
ElasticityLeft[k], ElasticityRight[k],WindowLeft[k], WindowRight[k],
MinLapseRate[k],MaxLapseRate[k] )
Price2[i] <- Price[k]*(1+x[i])*MarketPriceMean[i]
}
CumulRetentionCoef <- cumprod(1-R)
Result <-
ResultScenar(FirstYear,OptimYear,Price2,CumulRetentionCoef,Commission,IndividualExp,NbOfFutureSimYear,
IndivMeanLossAmount,
CPIRate,LossInflaRate,Alpha)
sum(Beta[1:NbOfFutureSimYear]*Result)
}
#########################################################
# Function calculating the optimal poistioning for a given problem year &
filtration :
#########################################################
YearFiltrationPbOptimPosition <- function(
PriorMarketPricePosition,MarketPricePhase,LRUpBound,LRLowBound,TransitCoefMean,TransitCoefStdDev,
LTMeanComputYearNb ,
FirstYear,LastYear,OptimYear,
CentralLapseRate, ElasticityLeft,WindowLeft, WindowRight, ElasticityRight,
Price,Commission,IndividualExp,IndivMeanLossAmount,
CPIRate,LossInflaRate,RiskFreeRate,RiskRate,Alpha)
{
Beta <- 1/cumprod(1+RiskFreeRate+RiskRate)
CumulLossInflaRate <- cumprod(1+LossInflaRate)
CumulCPIRate <- cumprod(1+CPIRate)
xmin <- Xmin(LapseType, ElasticityLeft, WindowLeft,
MinLapseRate,CentralLapseRate)
xmax <- Xmax(LapseType, ElasticityLeft, WindowLeft,
MaxLapseRate,CentralLapseRate)
kmin <- OptimYear - FirstYear + 1
kmax <- LastYear - FirstYear + 1
NbOfFutureSimYear <- LastYear - OptimYear + 1
PriorLR <- (LRUpBound*LRLowBound)^0.5
Phase <- 1
LRScenar <- LRScenar(PriorLR,Phase,LRUpBound
,LRLowBound,TransitCoefMean,TransitCoefStdDev,LTMeanComputYearNb)
MarketPriceLTMean <- mean( 1/LRScenar)
MarketPriceMatrix <-
MarketPriceMatrix(NbOfFutureSimYear+1,NbOfSimScenar,Phase,LRUpBound ,LRLowBound,
TransitCoefMean,TransitCoefStdDev,PriorMarketPricePosition,MarketPricePhase,MarketPriceLTMean)
MarketPriceMean <- apply(MarketPriceMatrix,2,mean)
Evaluate <- function(x)
{
ResultActuSum(FirstYear,OptimYear,x,NbOfFutureSimYear,
CentralLapseRate, ElasticityLeft,WindowLeft, WindowRight, ElasticityRight,
MarketPriceMean,Price,Commission,IndividualExp,IndivMeanLossAmount,
CPIRate,LossInflaRate,Beta,Alpha)
}
output <- genoud(Evaluate , nvars=NbOfFutureSimYear,Domains =
cbind(xmin[kmin : kmax],xmax[kmin : kmax]),
pop.size=5000,max=TRUE)
Solution <- rev(output$par)
MarketPriceMean <-rev(MarketPriceMean)
Solution
}
########################################################
# Problem parameters for a given year N & the coressponding filtration in year
N-1 :
########################################################
OptimYear <- 2017
# year of optimisation
PriorMarketPricePosition <- 1
#market price position in OptimYear-1
MarketPricePhase <- 1
#market price (excluding inflation) phase in OptimYear
#1: increasing, -1: decreasing
########################################################
Xoptim <- YearFiltrationPbOptimPosition(
PriorMarketPricePosition,MarketPricePhase,LRUpBound,LRLowBound,TransitCoefMean,TransitCoefStdDev,
LTMeanComputYearNb ,
FirstYear,LastYear,OptimYear,
CentralLapseRate, ElasticityLeft,WindowLeft, WindowRight, ElasticityRight,
Price,Commission,IndividualExp,IndivMeanLossAmount,
CPIRate,LossInflaRate,RiskFreeRate,RiskRate,Alpha)
Xoptim
Best Regards
thank you for your help
[[alternative HTML version deleted]]
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
