Régression contrainte dans R: coefficients positifs, somme à 1 et interception non nulle

8

J'ai le modèle que j'ai besoin d'estimer, avec et \ pi_k \ ge0 \ text {for} k \ geq 1 .

Y=π0+π1X1+π2X2+π3X3+ε,
∑kπk=1 for k≥1πk≥0 for k≥1

La réponse d' Elvis à une autre question résout ce problème dans le cas de π0=0 . Voici son code de cette solution:

   > library("quadprog");
   > X <- matrix(runif(300), ncol=3)
   > Y <- X %*% c(0.2,0.3,0.5) + rnorm(100, sd=0.2)
   > Rinv <- solve(chol(t(X) %*% X));
   > C <- cbind(rep(1,3), diag(3))
   > b <- c(1,rep(0,3))
   > d <- t(Y) %*% X  
   > solve.QP(Dmat = Rinv, factorized = TRUE, dvec = d, Amat = C, bvec = b, meq = 1)
   $solution
   [1] 0.2049587 0.3098867 0.4851546

   $value
   [1] -16.0402

   $unconstrained.solution
   [1] 0.2295507 0.3217405 0.5002459

   $iterations
   [1] 2 0

   $Lagrangian
   [1] 1.454517 0.000000 0.000000 0.000000

   $iact
   [1] 1

Comment puis-je ajuster ce code afin qu'il puisse estimer une interception?

Ceci a été mis en ligne ici parce que mon groupe dans ma mission est ennuyé de ne pas avoir encore estimé cette régression. Je répondrai à cette question ici si / quand les autres participants au forum y arriveront en premier.

r regression econometrics least-squares constrained-regression Communauté
la source

Réponses:

8

Vous avez juste besoin de jouer un peu avec les matrices impliquées. Ajoutez l'interception à X:

XX <- cbind(1,X)

Recalculez la Dmatrice utilisée dans solve.QP()(je préfère travailler directement avec cela pour éviter d'appeler solve():

Dmat <- t(XX)%*%XX

Recalculez davec le nouveau XX:

dd <- t(Y)%*%XX

Modifiez la matrice de contraintes en ajoutant une colonne zéro, car vous semblez ne pas avoir de contraintes sur l'ordonnée à l'origine (non?): 

Amat <- t(cbind(0,rbind(1,diag(3))))

Et enfin:

solve.QP(Dmat = Dmat, factorized = FALSE, dvec = dd, Amat = Amat, bvec = b, meq = 1)
Stephan Kolassa
la source
Merci Stephan. Je veux utiliser un bootstrap pour estimer les erreurs standard de ces estimations de coefficient. J'utiliserai une corrélation série et un bootstrap hétéro robuste. Avez-vous des opinions sur la légitimité de cette approche (le bootstrap en général, pas le type exact de bootstrap que je vais utiliser)?
En principe, il semble que le bootstrap soit légitime ici ... tant que vous gardez à l'esprit les contraintes sous lesquelles il a été calculé et communiquez clairement ces contraintes, mais vous alliez le faire, n'est-ce pas ;-)? Mais je suis loin d'être un expert en la matière; peut-être que quelqu'un de plus compétent peut commenter?
Stephan Kolassa
2
Je serais prudent à propos du bootstrap si l'un des , estimés est nul, car cette solution se trouve le long d'une frontière de l'espace des paramètres et les conditions de régularité utilisées pour justifier le bootstrap peuvent ne pas tenir. Mais c'est précisément le cas où cette solution est nécessaire: après tout, si tous les sont différents de zéro, alors on pourrait simplement utiliser OLS pour résoudre ce problème. πii≥1πi
whuber
@whuber: bon point. Bien sûr, l'autre objectif de cette chose est que les coefficients soient non seulement non négatifs, mais qu'ils totalisent également 1. Supposons que tous les coefficients soient bien éloignés de 0. Le bootstrap avec la contrainte de somme serait-il valide?
Stephan Kolassa
La somme limite les paramètres à un sous-collecteur de l'original (un sans limites ni coins). Ainsi, cela n'a pas le même effet que d'imposer une condition aux limites et ne devrait pas affecter l'estimation, l'amorçage, etc. Une autre façon de voir cela serait de réduire le nombre de paramètres d'une , par exemple en réécrivant et réécriture des conditions de non-négativité sous forme d'inégalités dans et uniquement. (Les conditions pour forcer toutes les valeurs à être inférieures ou égales à sont alors inutiles.)π3=1−(π1+π2)π1π21
whuber