2013-06-26 39 views
0
obj1<-function(monthly.savings, 
       success, 
       start.capital, 
       target.savings, 
       monthly.mean.return, 
       monthly.ret.std.dev, 
       monthly.inflation, 
       monthly.inf.std.dev, 
       n.obs, 
       n.sim=1000){ 

    req = matrix(start.capital, n.obs+1, n.sim) #matrix for storing target weight 

    monthly.invest.returns = matrix(0, n.obs, n.sim) 
    monthly.inflation.returns = matrix(0, n.obs, n.sim) 

    monthly.invest.returns[] = rnorm(n.obs * n.sim, mean = monthly.mean.return, sd = monthly.ret.std.dev) 
    monthly.inflation.returns[] = rnorm(n.obs * n.sim, mean = monthly.inflation, sd = monthly.inf.std.dev) 

    #for loop to be 
    for (a in 1:n.obs){ 
    req[a + 1, ] = req[a, ] * (1 + monthly.invest.returns[a,] - monthly.inflation.returns[a,]) + monthly.savings 
    } 

    ending.values=req[nrow(req),] 
    suc<-sum(ending.values>target.savings)/n.sim 
    value<-success-suc 

    return(abs(value)) 
} 

我有上面我想要最小化的目標函數。它試圖解決給定的成功概率所需的每月節約。考慮下面的輸入假設R優化給出的目標函數

success<-0.9 
start.capital<-1000000 
target.savings<-1749665 
monthly.savings=10000 
monthly.mean.return<-(5/100)/12 
monthly.ret.std.dev<-(3/100)/sqrt(12) 
monthly.inflation<-(5/100)/12 
monthly.inf.std.dev<-(1.5/100)/sqrt(12) 
monthly.withdrawals<-10000 
n.obs<-10*12 #years * 12 months in a year 
n.sim=1000 

我用下面的符號:

optimize(f=obj1, 
     success=success, 
     start.capital=start.capital, 
     target.savings=target.savings, 
     monthly.mean.return=monthly.mean.return, 
     monthly.ret.std.dev=monthly.ret.std.dev, 
     monthly.inflation=monthly.inflation, 
     monthly.inf.std.dev=monthly.inf.std.dev, 
     n.obs = n.obs, 
     n.sim = n.sim, 
     lower = 0, 
     upper = 10000, 
     tol = 0.000000001,maximum=F) 

我得到7875.03

因爲我從正態分佈取樣,輸出每一次會有所不同,但他們應該是相同的給予或採取幾個百分點。我遇到的問題是我無法任意指定上限。上面的例子的上限(10000)經過多次試驗後被挑選出櫻桃。如果說我把100000的上限(我知道是不合理的),它會返回這個數字,反對找到全球最低儲蓄。任何不正確的構思我的目標函數的想法?

感謝,

回答

4

你的函數並不總是返回相同的輸出給定輸入 很可能會提出幾個問題(它會創造出很多假的局部極小的)事實: 你能避免它們通過將函數(例如,set.seed(1)), 內的隨機數發生器 的種子或通過存儲隨機數並每次重複使用它們,或者通過使用低差異序列(例如,randtoolbox::sobol)來設置它們的種子。

由於它是一個變量的函數,因此您可以簡單地繪製它看看會發生什麼: 它在10,000之後有一個平臺 - 優化算法不能區分平臺和局部最優之間的 。

f <- function(x) { 
    set.seed(1) 
    obj1(x, 
     success    = success, 
     start.capital  = start.capital, 
     target.savings  = target.savings, 
     monthly.mean.return = monthly.mean.return, 
     monthly.ret.std.dev = monthly.ret.std.dev, 
     monthly.inflation = monthly.inflation, 
     monthly.inf.std.dev = monthly.inf.std.dev, 
     n.obs    = n.obs, 
     n.sim    = n.sim 
) 
} 
g <- Vectorize(f) 
curve(g(x), xlim=c(0, 20000)) 

您最初的問題其實不是一個最小化問題, 但求根問題,這是很容易。

obj2 <- function(monthly.savings) { 
    set.seed(1) 
    req = matrix(start.capital, n.obs+1, n.sim) 
    monthly.invest.returns <- matrix(0, n.obs, n.sim) 
    monthly.inflation.returns <- matrix(0, n.obs, n.sim) 
    monthly.invest.returns[] <- rnorm(n.obs * n.sim, mean = monthly.mean.return, sd = monthly.ret.std.dev) 
    monthly.inflation.returns[] <- rnorm(n.obs * n.sim, mean = monthly.inflation, sd = monthly.inf.std.dev) 
    for (a in 1:n.obs) 
    req[a + 1, ] <- req[a, ] * (1 + monthly.invest.returns[a,] - monthly.inflation.returns[a,]) + monthly.savings 
    ending.values <- req[nrow(req),] 
    suc <- sum(ending.values>target.savings)/n.sim 
    success - suc 
} 
uniroot(obj2, c(0, 1e6)) 
# [1] 7891.187