Как мне перевести эту проблему максимизации внутри уравнения в R? - PullRequest
Новая
       31

Как мне перевести эту проблему максимизации внутри уравнения в R?

0 голосов
/ 29 мая 2019

коллег-программистов. Я изучаю книгу о численных решениях для экономики (Judd 1998). Я пытаюсь воспроизвести проблему из той же книги в R, чтобы я мог использовать пакет optim, чтобы посмотреть, смогу ли я получить аналогичные результаты.

Проблема, установленная автором, заключается в следующем: one: , и его результаты были Эти .

Я попытался записать эту проблему в R, в результате чего получился фрагмент кода: 

DisutilityJudd <- function(L){
  if(L == 0){
    return(0)
  }else{
    return(0.1)
  }
}

AgentUtilityJudd <- function(w, L){
  (-exp(-2*w) + 1) - DisutilityJudd(L)
}

reservation.utility.judd <- AgentUtilityJudd(1, 1)

MaxEffortUtility <- function(w1, w2, L = 1){
  0.8 * AgentUtilityJudd(w1, L) + 0.2 * AgentUtilityJudd(w2, L)
}

LeastEffortUtility <- function(w1, w2, L = 0){
  0.4 * AgentUtilityJudd(w1, L) + 0.6 * AgentUtilityJudd(w2, L)
}

UtilityDifferenceJudd <- function(w1, w2){
  MaxEffortUtility(w1, w2) - LeastEffortUtility(w1, w2)
}

PenaltyFunctionJudd <- function(w1, w2, P = 100000){
  if(length(w1) == 2){
    y <- -1 *  (0.8 * (2 - w1[1]) - 0.2 * w1[2] - P * 
                  (pmax(0, -MaxEffortUtility(w1[1], w1[1]) - reservation.utility.judd))^2 -
                  P * (pmax(0, -UtilityDifferenceJudd(w1[1], w1[1])))^2)
  }else{
    y <- -1 * (0.8 * (2 - w1) - 0.2 * w2 - P * 
                 (pmax(0, -MaxEffortUtility(w1, w2) - reservation.utility.judd))^2 -
                 P * (pmax(0, -UtilityDifferenceJudd(w1, w2)))^2)
  }
  return(y)
}

Ошибок не было, но результаты, сгенерированные моим кодом, оказались далеко не такими, как я ожидал: 

optim(c(1.1, 0.5), PenaltyFunctionJudd)
$par
[1]  1.343909e+49 -2.370681e+51

$value
[1] -4.633849e+50

$counts
function gradient 
     501       NA 

$convergence
[1] 1

$message
NULL

Возможно, есть проблема с моей штрафной функцией. Я предполагаю, что это связано с функцией pmax. Может ли кто-нибудь помочь мне определить это? Спасибо, я ценю ваше внимание.

Редактировать: опечатка.

1 Ответ

0 голосов
/ 29 мая 2019

Я полагаю, вы имели в виду w1[2], когда if(length(w1) == 2) истинно.

Я изменил ваш код, не касаясь того, как вы определяете предыдущую функцию.Не ясно, если это ожидаемый результат: что означает IV (-1), это результат минус 1?сила, если 10? 

PenaltyFunctionJudd <- function(w1, w2, P = 1e5){

        if(length(w1) > 1){
                w2 <- w1[2]
                w1 <- w1[1]
        }
        # cat("length is 2 \n")

        y <- 0.8 * (2 - w1) - 0.2 * w2 - P *
                ( pmax(0, -MaxEffortUtility(w1, w2) - reservation.utility.judd) )^2 -
                P * ( pmax(0, -UtilityDifferenceJudd(w1, w2)) )^2

        # cat("pmax1 :", pmax(0, -MaxEffortUtility(w1, w2) - reservation.utility.judd), "\n")
        # cat("pmax2 :", pmax(0, -UtilityDifferenceJudd(w1, w2)), "\n")

        return(y)
}

optim(c(1.1, 0.5), PenaltyFunctionJudd, control = list(fnscale = -1) )
optim(c(11, 5), PenaltyFunctionJudd, method = "BFGS", control = list(fnscale = -1, maxit = 100) )

Вы можете использовать cat или print для проверки ваших значений (здесь я заметил, что некоторые Inf и 0 привели меня, чтобы заметить ошибку кода).

Дружественное предупреждение: при условии, что вы правильно определили предыдущую функцию, в оптимизации много нестабильности (проблема неправильно установлена? Требуется больше штрафов?).Действительно, при запуске дважды или более параметры алгоритма сильно колеблются ...

...