Comment puis-je ajuster une régression contrainte dans R afin que les coefficients totaux = 1?

36

Je vois une régression similaire similaire ici:

Régression linéaire contrainte à travers un point spécifié

mais mon exigence est légèrement différente. Il me faut les coefficients pour faire un total de 1. Plus précisément, je régresse les rendements d'une série de devises contre trois autres séries de devises, de sorte que les investisseurs puissent remplacer leur exposition à cette série par une combinaison de l'exposition aux 3 autres, mais leur les dépenses en espèces ne doivent pas changer, et de préférence (mais ce n'est pas obligatoire), les coefficients doivent être positifs.

J'ai essayé de rechercher une régression contrainte dans R et Google mais avec peu de chance.

Thomas Browne
la source
Êtes-vous sûr qu'il s'agit d'un problème de régression limité? En lisant la question, vous recherchez une relation de la forme (une série Forex) = (plus, je suppose, un quatrième terme représentant un taux de rendement sûr dominant). C'est indépendant de la décision d'investissement. Si un client veut investir capital dans y_4 en utilisant y_1 , y_2 et y_3 comme proxy, il investira simplement c \ beta_1 dans y_1 , c \ beta_2 dans y_2 et c \ beta_3 dans y_3y4β1y1+β2y2+β3y3cy 1 y 2 y 3 c β 1 y 1 c β 2y4y1y2y3cβ1y1cβ2 c β 3 y 3y2cβ3y3. Cela n'ajoute aucune complication particulière à la régression, n'est-ce pas?
whuber
C'est parce que si vous modélisez ceci, vous constaterez que B1 + B2 + B3> 1 dans de nombreux cas (ou <1 dans d'autres). En effet, la devise que l'on essaie de répliquer avec les descripteurs aura généralement une volatilité plus grande ou plus petite que les autres, de sorte que la régression vous donnera une pondération plus petite ou plus grande en réponse. Cela nécessite que l’investisseur soit ne soit pas pleinement investi, soit mette à profit ce que je ne souhaite pas. Quant au taux de rendement sûr no. Tout ce que nous essayons de faire est de répliquer series1 en utilisant d'autres variables. Étant un financier et non un statisticien, j'ai peut-être mal nommé ma question.
Thomas Browne
La raison pour inclure un terme pour un taux de rendement sûr est que parfois il aura un coefficient différent de zéro. Vraisemblablement, des instruments sécurisés (dépôts bancaires au jour le jour) sont accessibles à tous à faible coût, de sorte que quiconque l'ignorant en tant que composante de son panier d'investissement pourrait choisir des combinaisons non optimales. Maintenant, si les coefficients n’ajoutent pas à l’unité, alors quoi? Il suffit d’investir autant que vous le souhaitez dans les proportions estimées par la régression.
whuber
droit ..... aussi simple que ça. Merci. Je me sens un peu bête maintenant haha.
Thomas Browne
1
Pas idiot du tout. Le simple fait de poser cette question reflète un niveau élevé de réflexion. Je vérifiais simplement ma propre compréhension de votre question pour m'assurer que vous obteniez une réponse efficace. À votre santé.
whuber

Réponses:

35

Si je comprends bien, votre modèle est avec et . Vous devez minimiser sous réserve de ces contraintes. Ce type de problème est connu sous le nom de programmation quadratique .Σ k π k = 1 π k0 Σ i ( Y i - ( π 1 X i 1 + π 2 X i 2 + π 3 X i 3 ) ) 2

Y=π1X1+π2X2+π3X3+ε,
kπk=1πk0
i(Yi(π1Xi1+π2Xi2+π3Xi3))2

Voici quelques lignes de codes R donnant une solution possible ( sont les colonnes de , les valeurs vraies de sont 0.2, 0.3 et 0.5).π kX1,X2,X3Xπk

> 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

Je ne connais aucun résultat sur la distribution asymptotique des estimateurs, etc. Si quelqu'un a des indicateurs, je serais curieux d'en obtenir (si vous le souhaitez, je peux ouvrir une nouvelle question à ce sujet).

Elvis
la source
En fait question rapide. Ne devrais-je pas minimiser la variance plutôt que la somme? N’est-ce pas une régression qui minimise la variance du carré des erreurs?
Thomas Browne
6
C'est intelligent, Elvis, mais ne pourriez-vous pas accomplir la même chose simplement en reparamétrant la régression? Soit Cela équivaut à . Les estimations et les erreurs types de sont à calculer à partir des estimations et de la matrice var-covar de et . Y=α1X1+α2X2+(1α1α2)X3+εYX3=α1(X1X3)+α2(X2X3)+επiα1α2
whuber
6
@whuber Oui, mais avec des données plus bruitées ou avec certains des proches de , vous facilement la contrainte , qui constitue la partie "difficile" du problème. πk0πk>0
Elvis
2
Un coefficient positif vous indique d'acheter une devise étrangère; un coefficient négatif vous dit de le vendre. Si vous ne possédez pas déjà cette devise, vous devez l'emprunter pour la vendre ("vendre à découvert"). Parce que des emprunts illimités peuvent créer des problèmes, il existe des contraintes quant au montant de l'emprunt et à son mode de paiement ("exigences de marge", "coûts de constitution du capital" et "procédures d'évaluation à la valeur du marché"). Par conséquent, l’emprunt est possible mais est souvent évité sauf par les principaux acteurs des marchés ou à moins que cela ne confère d’importants avantages.
whuber
2
Un grand merci à tous pour toute l'aide. En fait, pour faire un commentaire sur les marchés des changes en général, il est plus facile de les vendre à découvert que les actions ou les obligations car il n’est pas nécessaire d’emprunter une action avant de vendre à découvert. On retourne simplement les monnaies dénominateur et numérateur. Ainsi, par exemple, vendre l'EURUSD et USDEUR sont des opérations tout à fait équivalentes du point de vue du département des risques, mais ce sont bien sûr des positions opposées. C'est pourquoi FX est un tel terrain de jeu pour les traders quant car ils n'ont pas à s'inquiéter des frictions directionnelles qui sont beaucoup plus importantes en actions
Thomas Browne
8

Comme mentionné par whuber, si vous êtes uniquement intéressé par les contraintes d'égalité, vous pouvez simplement utiliser la fonction standard lm () en réécrivant votre modèle:

Y=α+β1X1+β2X2+β3X3+ϵ=α+β1X1+β2X2+(1β1β2)X3+ϵ=α+β1(X1X3)+β2(X2X3)+X3+ϵ

Mais cela ne garantit pas que vos contraintes d'inégalité sont satisfaites! Cependant, dans ce cas, vous obtenez exactement le même résultat qu'avec l'exemple de programmation quadratique ci-dessus (en plaçant le X3 à gauche):

X <- matrix(runif(300), ncol=3)
Y <- X %*% c(0.2,0.3,0.5) + rnorm(100, sd=0.2)
X1 <- X[,1]; X2 <-X[,2]; X3 <- X[,3]
lm(Y-X3~-1+I(X1-X3)+I(X2-X3))
Matifou
la source
Dans le cas susmentionné de Matifou, comment éviter que le troisième coefficient soit négatif? Par exemple, si les coefficients optimaux pour et nous aurions cela ce qui implique ici que notre troisième coefficient est négatif et ne tient donc pas sur la base de notre régression souhaitée. β 2 = 0,5 ( 1 - β 1 - β 2 ) = - 0,25β1=0.75β2=0.5(1β1β2)=0.25
AS
1
Merci @AS pour l'avoir signalé. En effet, cette solution ne fonctionne que pour les contraintes d'égalité, pas pour celles d'inégalité. J'ai édité le texte en conséquence.
Matifou
1

Si je comprends bien votre modèle, vous cherchez à trouver tel que Σ[ ˉ b ]=1

x¯¯b¯=y¯
[b¯]=1

J'ai trouvé que le moyen le plus simple de traiter ce type de problèmes consiste à utiliser les propriétés associatives des matrices pour traiter en fonction d'autres variables.b¯

Par exemple, est une fonction de via le bloc de transformation . Dans votre cas, ci-dessous est . ici , nous pouvons séparer nos nowns et nknowns. ˉ c ¯ ¯ T c r1 ˉ b =[ k 0 k 1 k 2 ]= ¯ ¯ T c ˉ c =[ 1 0 0 0 1 0 - 1 - 1 1 ][ k 0 k 1 r ]kb¯c¯Tc¯¯r1

b¯=[k0k1k2]=Tc¯¯c¯=[100010111][k0k1r]
ku
c¯=[k0k1r]=Su¯¯cu¯+Sk¯¯ck¯=[100100][k0k1]+[001]r
Alors que je pouvais combiner les différentes transformations / blocs de séparation, cela devient lourd avec des modèles plus complexes. Ces blocs permettent de séparer les connus et les inconnus. ˉ ˉ v ¯ c u = ˉ w
x¯¯Tc¯¯(Su¯¯cu¯+Sk¯¯ck¯)=y¯v¯¯=x¯¯Tc¯¯Su¯¯w¯=y¯x¯¯Tc¯¯Sk¯¯ck¯
Enfin, le problème se présente sous une forme familière.
v¯¯cu¯=w¯
Augi Lynch
la source
1

Vieille question mais comme je suis confronté au même problème, j'ai pensé poster mon 2p ...

Utilisez la programmation quadratique comme suggéré par @Elvis mais en utilisant sqlincon à partir du paquet pracma . Je pense que l’avantage quadrpog::solve.QPest une interface utilisateur plus simple pour spécifier les contraintes. (En fait, lsqlinconc'est un wrapper autour solve.QP).

Exemple:

library(pracma)

set.seed(1234)

# Test data
X <- matrix(runif(300), ncol=3)
Y <- X %*% c(0.2, 0.3, 0.5) + rnorm(100, sd=0.2)

# Equality constraint: We want the sum of the coefficients to be 1.
# I.e. Aeq x == beq  
Aeq <- matrix(rep(1, ncol(X)), nrow= 1)
beq <- c(1)

# Lower and upper bounds of the parameters, i.e [0, 1]
lb <- rep(0, ncol(X))
ub <- rep(1, ncol(X))

# And solve:
lsqlincon(X, Y, Aeq= Aeq, beq= beq, lb= lb, ub= ub)

[1] 0.1583139 0.3304708 0.5112153

Mêmes résultats que ceux d'Elvis:

library(quadprog)
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

EDIT Pour essayer de répondre au commentaire de gung, voici une explication. sqlincon émule le lsqlin de matlab qui a une belle page d’aide. Voici les bits pertinents avec quelques modifications (mineures):

XMatrice multiplicateur, spécifiée comme une matrice de doubles. C représente le multiplicateur de la solution x dans l'expression C * x - Y. C est M-by-N, où M est le nombre d'équations et N le nombre d'éléments de x.

YVecteur constant, spécifié en tant que vecteur de doubles. Y représente le terme constant additif dans l'expression C * x - Y. Y est M par 1, où M est le nombre d'équations.

Aeq: Matrice de contrainte d'égalité linéaire, spécifiée comme une matrice de doubles. Aeq représente les coefficients linéaires dans les contraintes Aeq * x = beq. Aeq a la taille Meq-by-N, où Meq est le nombre de contraintes et N le nombre d'éléments de x

beqVecteur de contrainte d'égalité linéaire, spécifié en tant que vecteur de doubles. beq représente le vecteur constant dans les contraintes Aeq * x = beq. beq a la longueur Meq, où Aeq est Meq-by-N.

lbLimites inférieures, spécifiées en tant que vecteur de doubles. lb représente la limite inférieure élément par élément dans lb ≤ x ≤ ub.

ubLimites supérieures, spécifiées en tant que vecteur de doubles. ub représente la limite supérieure par élément dans lb ≤ x ≤ ub.

dariober
la source