[R] 12th Root of a Square (Transition) Matrix

David Winsemius dwinsemius at comcast.net
Fri Nov 4 22:37:50 CET 2011 

On Nov 4, 2011, at 4:04 PM, Peter Langfelder wrote:

> Is it just me or are you confusing the 12th root of a matrix with
> taking the 12th root of each entry?

I think I got confused as well. Thanks for clarifying.

> Because your formula involving the
> eigenvectors and eigenvalues calculates the 12th root of the matrix,
> while
>
> round(nth_root^12,4)
>
> will print out a matrix whose components are powers of 12 of the
> nth_root, which is very different. To find the 12th power of a matrix,
> you can either search for an appropriate function, or do
>
> res = nth_root;
> for (p in 1:11)
>  res = res %*% nth_root


The 12th (matrix) root of M: e^( 1/n * log(M) )

 > require(Matrix)
 > M1.12 <- expm( (1/12)*logm(M) )

 > res = M1.12; nth_root=M1.12
 > for (p in 1:11)
+  res = res %*% nth_root

# Check accuracy
 > round(res, 4)
        [,1]   [,2]   [,3]   [,4]   [,5]   [,6]   [,7]   [,8]
[1,] 0.9081 0.0833 0.0068 0.0006 0.0008 0.0002 0.0001 0.0001
[2,] 0.0070 0.9065 0.0779 0.0064 0.0006 0.0013 0.0002 0.0001
[3,] 0.0009 0.0227 0.9105 0.0552 0.0074 0.0026 0.0001 0.0006
[4,] 0.0002 0.0033 0.0595 0.8593 0.0530 0.0117 0.0112 0.0018
[5,] 0.0003 0.0014 0.0067 0.0773 0.8053 0.0884 0.0100 0.0106
[6,] 0.0001 0.0011 0.0024 0.0043 0.0648 0.8346 0.0407 0.0520
[7,] 0.0021 0.0000 0.0022 0.0130 0.0238 0.1124 0.6486 0.1979
[8,] 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 1.0000

 > M
        AAA     AA      A    BBB     BB      B    CCC      D
AAA 0.9081 0.0833 0.0068 0.0006 0.0008 0.0002 0.0001 0.0001
AA  0.0070 0.9065 0.0779 0.0064 0.0006 0.0013 0.0002 0.0001
A   0.0009 0.0227 0.9105 0.0552 0.0074 0.0026 0.0001 0.0006
BBB 0.0002 0.0033 0.0595 0.8593 0.0530 0.0117 0.0112 0.0018
BB  0.0003 0.0014 0.0067 0.0773 0.8053 0.0884 0.0100 0.0106
B   0.0001 0.0011 0.0024 0.0043 0.0648 0.8346 0.0407 0.0520
CCC 0.0021 0.0000 0.0022 0.0130 0.0238 0.1124 0.6486 0.1979
D   0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 1.0000

Rather good agreement to 4 decimal places anyway.


>
> then compare res to your original matrix M.
>
> Peter
>
> On Fri, Nov 4, 2011 at 3:34 AM, clangkamp <christian.langkamp at gmxpro.de 
> > wrote:
>> I have tried this method, but the result is not working, at least  
>> not as I
>> expect:
>> I used the CreditMetrics package transition matrix
>> 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)
>>
>> then followed through with the steps:
>>
>> nth_root <- X %*% L_star %*% X_inv
>>
>> But the check (going back 12 to the power again) doesn't yield the  
>> original
>> matrix. Now some rounding errors can be expected, but I didn't  
>> expect a
>> perfectly diagonal matrix, when the initial matrix isn't diagonal  
>> at all.
>>> round(nth_root^12,4)
>>       [,1]   [,2]   [,3]   [,4]   [,5]   [,6]   [,7] [,8]
>> [1,] 0.9078 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000    0
>> [2,] 0.0000 0.9053 0.0000 0.0000 0.0000 0.0000 0.0000    0
>> [3,] 0.0000 0.0000 0.9079 0.0000 0.0000 0.0000 0.0000    0
>> [4,] 0.0000 0.0000 0.0000 0.8553 0.0000 0.0000 0.0000    0
>> [5,] 0.0000 0.0000 0.0000 0.0000 0.7998 0.0000 0.0000    0
>> [6,] 0.0000 0.0000 0.0000 0.0000 0.0000 0.8285 0.0000    0
>> [7,] 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.6457    0
>> [8,] 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000    1
>>
>> Any takers
>>
>>
>> -----
>> Christian Langkamp
>> christian.langkamp-at-gmxpro.de
>>
>> --
>> View this message in context: http://r.789695.n4.nabble.com/12th-Root-of-a-Square-Transition-Matrix-tp2259736p3989618.html
>> Sent from the R help mailing list archive at Nabble.com.
>>
>> ______________________________________________
>> R-help at r-project.org 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.
>>
>
>
>
> -- 
> Sent from my Linux computer. Way better than iPad :)
>
> ______________________________________________
> R-help at r-project.org 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.

David Winsemius, MD
Heritage Laboratories
West Hartford, CT

More information about the R-help mailing list