Come impostare la somma dei parametri su 1 nell'ottimizzazione vincolata

Question

Ecco il codice (mi dispiace se è così lungo, ma è stato il primo esempio che ho avuto); Sto utilizzando l'esempio CVaR da CreditMetrics pacchetto di A. Wittmann e DEoptim risolutore di ottimizzare:

library(CreditMetrics) 
library(DEoptim) 

N <- 3 
n <- 100000 
r <- 0.003 
ead <- rep(1/N,N) 
rc <- c("AAA", "AA", "A", "BBB", "BB", "B", "CCC", "D") 
lgd <- 0.99 
rating <- c("BBB", "AA", "B") 
firmnames <- c("firm 1", "firm 2", "firm 3") 
alpha <- 0.99 

# correlation matrix 
rho <- matrix(c( 1, 0.4, 0.6, 
        0.4, 1, 0.5, 
        0.6, 0.5, 1), 3, 3, dimnames = list(firmnames, firmnames), 
       byrow = TRUE) 

# one year empirical migration matrix from standard&poors website 
rc <- c("AAA", "AA", "A", "BBB", "BB", "B", "CCC", "D") 
M <- matrix(c(90.81, 8.33, 0.68, 0.06, 0.08, 0.02, 0.01, 0.01, 
       0.70, 90.65, 7.79, 0.64, 0.06, 0.13, 0.02, 0.01, 
       0.09, 2.27, 91.05, 5.52, 0.74, 0.26, 0.01, 0.06, 
       0.02, 0.33, 5.95, 85.93, 5.30, 1.17, 1.12, 0.18, 
       0.03, 0.14, 0.67, 7.73, 80.53, 8.84, 1.00, 1.06, 
       0.01, 0.11, 0.24, 0.43, 6.48, 83.46, 4.07, 5.20, 
       0.21,  0, 0.22, 1.30, 2.38, 11.24, 64.86, 19.79, 
       0,  0,  0,  0,  0,  0,  0, 100 
)/100, 8, 8, dimnames = list(rc, rc), byrow = TRUE) 

cm.CVaR(M, lgd, ead, N, n, r, rho, alpha, rating) 

y <- cm.cs(M, lgd)[which(names(cm.cs(M, lgd)) == rating)] 

Ora scrivo la mia funzione ...

fun <- function(w) { 
    # ... 
    - (t(w) %*% y - r)/cm.CVaR(M, lgd, ead = w, N, n, r, 
          rho, alpha, rating) 
} 

... e voglio ottimizzare iT:

DEoptim(fn = fun, lower = rep(0, N), upper = rep(1, N), 
     control = DEoptim.control()) 

Potete dirmi cosa devo inserire in # ... per fare sum(w) = 1 durante l'ottimizzazione?

Qui di seguito vi mostro i risultati di ottimizzazione in base alle punte di flodel:

# The first trick is to include B as large number to force the algorithm to put sum(w) = 1 

fun <- function(w) { 
    - (t(w) %*% y - r)/cm.CVaR(M, lgd, ead = w, N, n, r, rho, alpha, rating) + 
    abs(10000 * (sum(w) - 1)) 
} 

DEoptim(fn = fun, lower = rep(0, N), upper = rep(1, N), 
     control = DEoptim.control()) 

$optim$bestval 
[1] -0.05326055 

$optim$bestmem 
par1  par2  par3 
0.005046258 0.000201286 0.994752456 

parsB <- c(0.005046258, 0.000201286, 0.994752456) 

> fun(parsB) 
      [,1] 
[1,] -0.05326089 

... e ...

Come si può vedere, il primo trucco funziona meglio in che trova un risultato che è più piccolo del secondo Sfortunatamente sembra che impieghi più tempo.

# The second trick needs you use w <- w/sum(w) in the function itself 

fun <- function(w) { 
    w <- w/sum(w) 
    - (t(w) %*% y - r)/cm.CVaR(M, lgd, ead = w, N, n, r, rho, alpha, rating) #+ 
    #abs(10000 * (sum(w) - 1)) 
} 

DEoptim(fn = fun, lower = rep(0, N), upper = rep(1, N), 
     control = DEoptim.control()) 

$optim$bestval 
[1] -0.0532794 

$optim$bestmem 
par1   par2   par3 
1.306302e-15 2.586823e-15 9.307001e-01 

parsC <- c(1.306302e-15, 2.586823e-15, 9.307001e-01) 
parC <- parsC/sum(parsC) 

> fun(parC) 
      [,1] 
[1,] -0.0532794 

Qualsiasi commento?

Devo aumentare il numero di iterazioni a causa di una funzione "troppo stocastica" da ottimizzare?

Answer 1

Prova:

w <- w/sum(w) 

e se DEoptim ti dà una soluzione ottimale w* tale che sum(w*) != 1 poi w*/sum(w*) dovrebbe essere la vostra soluzione ottimale.

Un altro approccio è quello di risolvere tutte le variabili tranne una. Sappiamo il valore dell'ultimo variabile deve essere 1 - sum(w) così nel corpo della funzione, avere:

w <- c(w, 1-sum(w)) 

e fare lo stesso per la soluzione ottimale restituito da DEoptim: w* <- c(w*, 1-sum(w*))

Entrambe le soluzioni richiedono che si riformulare il problema in un'ottimizzazione non vincolata (senza contare per i limiti delle variabili) in modo da poter utilizzare DEoptim; che ti costringe a fare un piccolo lavoro extra al di fuori di DEoptim per recuperare la soluzione al problema originale.

In risposta al tuo commento, se vuoi che DEoptim ti dia la risposta corretta subito (cioè senza la necessità di una post-trasformazione), potresti anche provare ad includere un costo di penalità per la tua funzione obiettivo: ad esempio aggiungi B * abs(sum(w)-1) dove B è un numero arbitrario di grandi dimensioni, quindi sum(w) sarà forzato a 1.

Answer 2

Penso che dovresti aggiungere una penalità per ogni deviazione da uno. Aggiungi al tuo problema di minimizzazione il termine +(sum(weights) - 1)^2 * 1e10. Dovresti vedere che questa enorme penalità costringerà i pesi a sommarsi a 1!