R help - Dec 2004 - Re: Help : generating correlation matrix with a particular

************************************************************
Important: We would  draw your  attention to the  notices at 
the bottom of this  e-mail, particularly  before opening and 
reviewing any file attachment(s).
************************************************************

Here is some code we have used.

a<-array(c(1,.9,.7,.9,1,.3,.7,.3,1),dim=c(3,3))
a
s<-eigen(a)$vectors
l<-diag(eigen(a)$values)
l[l<0]<-0
b<-s%*%sqrt(l)
for(i in 1:nrow(b)){b[i,]<-b[i,]/sqrt(sum(b[i,]^2))}
ap<-b%*%t(b)
ap


It is based on a paper by Rebonato etal that formerly was at
www.rebonato.com/correlationmatrix.pdf.
Unfortunately the website has disappeared.

Best regards,
Herb


Herbert G. Desson, ACAS, MAAA

Actuary
JLT Risk Solutions
6 Crutched Friars
London EC3N 2PH

phone:  +44 (0)20 7528 4702
fax:       +44 (0)20 7558 3785



>Message: 2
>Date: 12 Dec 2004 14:58:38 +0100
>From: Peter Dalgaard <p.dalgaard at biostat.ku.dk>
>Subject: Re: [R] Help : generating correlation matrix with a
>	particular	structure
>To: Siew Leng TENG <siewlengteng at yahoo.com>
>Cc: r-help at stat.math.ethz.ch
>Message-ID: <x2sm6bd3e9.fsf at biostat.ku.dk>
>Content-Type: text/plain; charset=us-ascii
>
>Siew Leng TENG <siewlengteng at yahoo.com> writes:
>
>> Hi,
>> 
>> I would like to generate a correlation matrix with a
>> particular structure. For example, a 3n x 3n matrix :
>> A_(nxn)   aI_(nxn)  bI_(nxn)
>> aI_(nxn)  A_(nxn)   cI_(nxn)
>> aI_(nxn)  cI_(nxn)  A_(nxn)
>> 
>> where
>> - A_(nxn) is a *specified* symmetric, positive
>> definite nxn matrix.
>> - I_(nxn) is an identity matrix of order n
>> - a, b, c are (any) real numbers
>> 
>> Many attempts have been unsuccessful because a
>> resulting matrix with any a, b, c may not be a
>> positive definite one, and hence cannot qualify as a
>> correlation matrix. Trying to first generate a
>> covariance matrix however, does not guarantee a
>> corresponding correlation matrix with the above
>> structure.
>
>Er, a correlation matrix *is* a covariance matrix with 1 down the
>diagonal... 
>
>You need to sort out the parametrization issues. What you're trying to
>achieve is quite hard. Consider the simpler case of two blocks and
>n=2; what you're asking for is a covariance matrix of the form
>
>1 r a 0
>r 1 0 a
>a 0 1 r
>0 a r 1
>
>so if this is the correlation matrix of (X1,Y1,X2,Y2) you want
>
>X1 and Y1 correlated 
>X2 and Y2 correlated
>X1 and X2 correlated
>Y1 and Y2 correlated
>
>but
>
>X1 and Y2 uncorrelated
>Y1 and X2 uncorrelated
>
>
>One approach is to work out the conditional variance of (X2,Y2) given
>(X1,Y1) and check for positive semidefiniteness. You do the math...
>
>(Some preliminary experiments suggest that the criterion could be
>abs(a)+abs(r) <= 1, but don't take my word for it)
>
>> R-version used :
>> ---------------
>> Windows version
>> R-1.8.1
>> Running on Windows XP
>
>You might want to upgrade, but it might not do anything for you in
>this respect.
>
>-- 
>   O__  ---- Peter Dalgaard             Blegdamsvej 3  
>  c/ /'_ --- Dept. of Biostatistics     2200 Cph. N   
> (*) \(*) -- University of Copenhagen   Denmark      Ph: (+45) 35327918
>~~~~~~~~~~ - (p.dalgaard at biostat.ku.dk)             FAX: (+45) 35327907


************************************************************
JLT Risk Solutions Ltd
6 Crutched Friars, London EC3N 2PH. Co Reg No 1536540
Tel: (44) (0)20 7528 4000   Fax: (44) (0)20 7528 4500
http://www.jltgroup.com
Lloyd's Broker.  Regulated by the General Insurance
Standards Council
------------------------------------------------------------
The content of this e-mail (including any attachments) as 
received may not be the same as sent. If you consider that 
the content is material to the formation or performance of 
a contract or you are otherwise relying upon its accuracy, 
you should consider requesting a copy be sent by facsimile 
or normal mail.  In any event, please check this message
and any identified file attachment(s) upon receipt and notify
the sender immediately if there is any manifest transmission
error, omission or corruption. This does not change or reduce
any party's duty of utmost good faith when contracting for
insurance or reinsurance. The information in this e-mail is 
confidential and may be legally privileged. If you are not 
the intended recipient, please notify the sender immediately 
and then delete this e-mail entirely - you must not retain, 
copy, distribute or use this e-mail for any purpose or 
disclose any of its content to others.

Opinions, conclusions and other information in this e-mail 
that do not relate to the official business of JLT Risk 
Solutions Ltd shall be understood as neither given nor 
endorsed by it.  Please note we intercept and monitor 
incoming / outgoing e-mail and therefore you should neither 
expect nor intend any e-mail to be private in nature.

We have checked this e-mail for viruses and other harmful 
components and believe but not guarantee it virus-free prior 
to leaving our computer system.  However, you should satisfy 
yourself that it is free from harmful components, as we do 
not accept responsibility for any loss or damage it may 
cause to your computer systems.

>>>>> "Herbert" == Herbert Desson <Herbert_Desson at
jltgroup.com>
>>>>>     on Mon, 13 Dec 2004 12:24:10 -0000 writes:

    Herbert> Here is some code we have used.

    Herbert> a<-array(c(1,.9,.7,.9,1,.3,.7,.3,1),dim=c(3,3))
    Herbert> a
    Herbert> s<-eigen(a)$vectors
    Herbert> l<-diag(eigen(a)$values)
    Herbert> l[l<0]<-0
    Herbert> b<-s%*%sqrt(l)
    Herbert> for(i in 1:nrow(b)){b[i,]<-b[i,]/sqrt(sum(b[i,]^2))}
    Herbert> ap<-b%*%t(b)
    Herbert> ap

This code does the same thing as my (simplistic, but slightly more
general) function posdefify()  in package "sfsmisc" :

  a <- matrix(c(1,.9,.7,.9,1,.3,.7,.3,1), 3)

  install.packages("sfsmisc")
  library(sfsmisc)

  posdefify(a)

gives
          [,1]      [,2]      [,3]
[1,] 1.0000000 0.8940242 0.6963190
[2,] 0.8940242 1.0000000 0.3009691
[3,] 0.6963190 0.3009691 1.0000000


    Herbert> It is based on a paper by Rebonato etal that formerly was at
    Herbert> www.rebonato.com/correlationmatrix.pdf.
    Herbert> Unfortunately the website has disappeared.

The idea is very simple and has been re-invented many times as
far as I know.

More sophisticated methods for "posdefiying" a matrix exist in
other places. Given symmetrix matrix A, they try to find the
matrix Ap, positive definite, such  ||A - Ap||  is minimal.
The eigen-value based simple solution that you've used above
and I've also coded in posdefify(), is not the same one would
get for `usual' matrix norms  || . || 

[[NB:  posdefify() also has a 2nd method the implementation of
       which has an embarassing bug. The next version of
       sfsmisc, due in a day or two, will have it fixed.
]]

Does anyone know of rigorous mathematical results in this
regard?

Martin Maechler, ETH Zurich