thr3ads.net - R help - [R] The quadprog package [Sep 2007]

If this information is useful, please help other people find it:
Share via:

thomas.schwander at mvv.de

2007-Sep-03 08:08 UTC

[R] The quadprog package

Hi everybody,

I'm using Windows XP Prof, R 2.5.1 and a Pentium 4 Processor.
Now, I want to solve a quadratic optimization program (Portfolio Selection) with
the quadprog package

I want to minimize (\omega'%*%\Sigma%*%\omega)
Subject to
	(1) \iota' %*% \omega = 1 (full investment)
	(2) R'%*%\omega = \mu (predefined expectation value)
	(3) \omega \ge 0 (no short sales).

Where 
	\omega is a Nx1 vector of the weights invested in each asset
	\Sigma is the NxN variance-covariance matrix
	\iota is a 1xN vector of 1's
	R' is a Nx1 vector of the expected returns
	and
	\mu is a number, the postualted return of the investor

I've done the following code but it doesn't make what I want it to do.
The weights aren't all positive and the \mu isn't reached. What's
wrong with my code?

Require(quadprog)
Dmat<-diag(1,7,7)
# muss als quadratische Matrix eingegeben werden
Dmat
dvec<-matrix(0,7,1) # muss als Spaltenvektor eingegeben werden
dvec
mu<-0 # (in Mio. €)
bvec<-c(1,mu,matrix(0,7,1)) # muss als Spaltenvektor eingegeben werden
bvec
mu_r<-c(19.7,33.0,0.0,49.7, 82.5, 39.0,11.8)	
Amat<-matrix(c(matrix(1,1,7),7*mu_r,diag(1,7,7)),9,7,byrow=T) 
# muss als Matrix angegeben werden, wie sie wirklich ist
Amat
meq<-2
loesung<-solve.QP(Dmat,dvec,Amat=t(Amat),bvec=bvec,meq=2)
loesung

# Überprüfen, ob System richtig gelöst wurde
loesung$solution %*% mu_r
sum(loesung$solution)
for (i in 1:7){
	a<-loesung$solution[i]>=0
	print(a)
}

Thanks in advance for your answers.

______________________________

Thomas Schwander

MVV Energie
Konzern-Risikocontrolling

Telefon 0621 - 290-3115
Telefax 0621 - 290-3664

E-Mail: Thomas.Schwander@mvv.de .  Internet: www.mvv.de
MVV Energie AG . Luisenring 49 . 68159 Mannheim
Handelsregister-Nr. HRB 1780, Amtsgericht Mannheim
Vorsitzender des Aufsichtsrates: Herr Dr. Peter Kurz
Vorstand: Dr. Rudolf Schulten (Vorsitzender) .  Matthias Brückmann . Dr. Werner
Dub . Hans-Jürgen Farrenkopf




	[[alternative HTML version deleted]]

Berwin A Turlach

2007-Sep-03 09:15 UTC

head link

[R] The quadprog package

G'day Thomas,

On Mon, 3 Sep 2007 10:08:08 +0200
<thomas.schwander at mvv.de> wrote:
> What's wrong with my code?
> 
> Require(quadprog)
library(quadprog) ? :)
> Dmat<-diag(1,7,7)
> # muss als quadratische Matrix eingegeben werden
> Dmat
> dvec<-matrix(0,7,1) # muss als Spaltenvektor eingegeben werden
> dvec
> mu<-0 # (in Mio. ?)
> bvec<-c(1,mu,matrix(0,7,1)) # muss als Spaltenvektor eingegeben werden
> bvec
> mu_r<-c(19.7,33.0,0.0,49.7, 82.5, 39.0,11.8)	
> Amat<-matrix(c(matrix(1,1,7),7*mu_r,diag(1,7,7)),9,7,byrow=T) 
> # muss als Matrix angegeben werden, wie sie wirklich ist
> Amat
> meq<-2
> loesung<-solve.QP(Dmat,dvec,Amat=t(Amat),bvec=bvec,meq=2)
With the commands above, I get on my system:
> loesung<-solve.QP(Dmat,dvec,Amat=t(Amat),bvec=bvec,meq=2)Error in solve.QP(Dmat, dvec, Amat = t(Amat), bvec = bvec, meq = 2) : 
	constraints are inconsistent, no solution!

Which is a bit strange, since Dmat is the identity matrix, so there
should be little room for numerical problems.

OTOH, it is known that the Goldfarb-Idnani algorithm can have problem
in rare occasions if the problem is "ill-scaled"; see the work of
Powell. 

Changing the definition of Amat to
> Amat<-matrix(c(matrix(1,1,7),mu_r,diag(1,7,7)),9,7,byrow=T) 
leads to a successful call to solve.QP.  And since mu=0, I really
wonder why you scaled up the vector mu_r by a factor of 7 when putting
it into Amat.... :)
> loesung<-solve.QP(Dmat,dvec,Amat=t(Amat),bvec=bvec,meq=2)
Now, with the solution I obtained the rest of your commands show:
> loesung$solution %*% mu_r             [,1]
[1,] 6.068172e-15> sum(loesung$solution)
[1] 1> for (i in 1:7){	a<-loesung$solution[i]>=0
	print(a)
}
[1] TRUE
[1] TRUE
[1] TRUE
[1] TRUE
[1] FALSE
[1] FALSE
[1] TRUE
 
But:
>  for (i in 1:7){	a<-loesung$solution[i]>=- .Machine$double.eps*1000
	print(a)
}

[1] TRUE
[1] TRUE
[1] TRUE
[1] TRUE
[1] TRUE
[1] TRUE
[1] TRUE
>  for (i in 1:7) print(loesung$solution[i])[1] 0
[1] 0
[1] 1
[1] 0
[1] -4.669169e-17
[1] -1.436586e-17
[1] 8.881784e-16

This is a consequence of finite precision arithmetic.  Just as you
should not compare to numeric numbers directly for equality but rather
that the absolute value of their difference is smaller than an
appropriate chosen threshold, you should not check whether a number is
bigger or equal to zero, but whether it is bigger or equal to an
appropriately chosen negative threshold.  More information are given in
FAQ 7.31.

HTH.

Cheers,

	Berwin

=========================== Full address ============================Berwin A
Turlach                            Tel.: +65 6515 4416 (secr)
Dept of Statistics and Applied Probability        +65 6515 6650 (self)
Faculty of Science                          FAX : +65 6872 3919       
National University of Singapore     
6 Science Drive 2, Blk S16, Level 7          e-mail: statba at nus.edu.sg
Singapore 117546                    http://www.stat.nus.edu.sg/~statba

Moshe Olshansky

2007-Sep-04 00:26 UTC

head link

[R] The quadprog package

Hi Thomas,

On my computer the solution is (0,0,1,0,0,0,0) (within
machine accuracy) and it satisfies the constraints.

--- thomas.schwander at mvv.de wrote:
> Hi everybody,
> 
> I'm using Windows XP Prof, R 2.5.1 and a Pentium 4
> Processor.
> Now, I want to solve a quadratic optimization
> program (Portfolio Selection) with the quadprog
> package
> 
> I want to minimize (\omega'%*%\Sigma%*%\omega)
> Subject to
> 	(1) \iota' %*% \omega = 1 (full investment)
> 	(2) R'%*%\omega = \mu (predefined expectation
> value)
> 	(3) \omega \ge 0 (no short sales).
> 
> Where 
> 	\omega is a Nx1 vector of the weights invested in
> each asset
> 	\Sigma is the NxN variance-covariance matrix
> 	\iota is a 1xN vector of 1's
> 	R' is a Nx1 vector of the expected returns
> 	and
> 	\mu is a number, the postualted return of the
> investor
> 
> I've done the following code but it doesn't make
> what I want it to do. The weights aren't all
> positive and the \mu isn't reached. What's wrong
> with my code?
> 
> Require(quadprog)
> Dmat<-diag(1,7,7)
> # muss als quadratische Matrix eingegeben werden
> Dmat
> dvec<-matrix(0,7,1) # muss als Spaltenvektor
> eingegeben werden
> dvec
> mu<-0 # (in Mio. ???)
> bvec<-c(1,mu,matrix(0,7,1)) # muss als Spaltenvektor
> eingegeben werden
> bvec
> mu_r<-c(19.7,33.0,0.0,49.7, 82.5, 39.0,11.8)	
>
Amat<-matrix(c(matrix(1,1,7),7*mu_r,diag(1,7,7)),9,7,byrow=T)> 
> # muss als Matrix angegeben werden, wie sie wirklich
> ist
> Amat
> meq<-2
>
loesung<-solve.QP(Dmat,dvec,Amat=t(Amat),bvec=bvec,meq=2)> loesung
> 
> # ??berpr??fen, ob System richtig gel??st wurde
> loesung$solution %*% mu_r
> sum(loesung$solution)
> for (i in 1:7){
> 	a<-loesung$solution[i]>=0
> 	print(a)
> }
> 
> Thanks in advance for your answers.
> 
> ______________________________
> 
> Thomas Schwander
> 
> MVV Energie
> Konzern-Risikocontrolling
> 
> Telefon 0621 - 290-3115
> Telefax 0621 - 290-3664
> 
> E-Mail: Thomas.Schwander at mvv.de .  Internet:
> www.mvv.de
> MVV Energie AG . Luisenring 49 . 68159 Mannheim
> Handelsregister-Nr. HRB 1780, Amtsgericht Mannheim
> Vorsitzender des Aufsichtsrates: Herr Dr. Peter Kurz
> Vorstand: Dr. Rudolf Schulten (Vorsitzender) . 
> Matthias Br??ckmann . Dr. Werner Dub . Hans-J??rgen
> Farrenkopf 
> 
> 
> 
> 
> 	[[alternative HTML version deleted]]
> 
> > ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained,
> reproducible code.
>

Reasonably Related Threads

Search for more possibly parallel threads

R help - Sep 2007 - The quadprog package

[R] The quadprog package

[R] The quadprog package

[R] The quadprog package

Reasonably Related Threads