Thank you SARAH
________________________________
De : Sarah Goslee <sarah.gos...@gmail.com>
rg>
Envoyé le : Vendredi, 31 Juillet 2009, 15h21mn 50s
Objet : Re: [R] what meaning missing value True /False needed
You are far more likely to get helpful volunteers to answer your
questions if you
distill your problem down to the smallest point. That's a lot of code to wade
through. But since I was waiting for something else to run, I took a stab at it.
It also would have been helpful if you'd given the *complete* error message,
since it tells you where exactly the problem is.
In this case, the problem is in LapseRate(), where PriceFactor is compared to
WindowLeft (as the actual error message says). I added two lines using cat()
to print those two variables each time LapseRate() is called, and here's what
I saw:
Maximization Problem.
PriceFactor: 0.9352139
WindowLeft: 1
PriceFactor: 1.627354
WindowLeft: 1
PriceFactor: NA
WindowLeft: 1
Error in if (PriceFactor < WindowLeft) r <- CentralLapseRate * (1 +
ElasticityLeft/WindowLeft * :
missing value where TRUE/FALSE needed
So indeed, LapseRate() is comparing a NA value to something else,
resulting in NA,
just as the error message stated.
Figuring out why PriceFactor is NA after two iterations is left as an
exercise for the
reader.
Sarah
rote:
> 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]]
>
>
--
Sarah Goslee
http://www.functionaldiversity.org
[[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.
