From: Duncan Murdoch <murdoch.duncan_at_gmail.com>

Date: Fri, 18 Jun 2010 06:13:23 -0400

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. Received on Fri 18 Jun 2010 - 10:21:47 GMT

Date: Fri, 18 Jun 2010 06:13:23 -0400

On 18/06/2010 2:01 AM, Corey Gallon wrote:

> Dear R-tisans,

*>
**> I am trying to calculate the 12th root of a transition (square) matrix, but can't seem to obtain an accurate result. I realize that this post is laced with intimations of quantitative finance, but the question is both R-related and broadly mathematical. That said, I'm happy to post this to R-SIG-Finance if I've erred in posting this to the general list.
**>
**> I've pulled down an annual transition matrix from the latest Moody's Corporate Default Study, and I'm using this (with the default row added manually) as the basis for this calculation. (I've pasted the dput of the resulting matrix below.) According to Hull, Appendix E [1], an arbitrary root of a square matrix (A) can be calculated by multiplying the inverse matrix of eigenvectors (X-inv) by the nth-root of diagonalized matrix of eigenvalues (Lambda-star) by the matrix of eigenvectors (X) -- all of these eigenvectors(values) being calculated from the matrix for which one wishes to calculate the nth root. The equation is as follows:
**>
**> A = X-inv %*% Lambda-star %*% X
**>
*

This is wrong: you've swapped X and X-inv. So your final line below should be

nth_root <- X %*% L_star %*% X_inv

Duncan Murdoch

> I've written the code below to implement this, but the result doesn't seem to be correct. (I can't raise the resulting matrix to the 12th power to calculate the original matrix.) I believe that the reason for this is the order in which R returns the eigenvalues (i.e. a vector in descending order) and the order in which I've created the matrix of eigenvectors, but I may be wrong in this suspicion.

*>
**> I defer to the collective wisdom of the community, and hope that minds greater than mine may provide insight.
**>
**>
**>
**> Cheers,
**>
**> Corey
**>
**>
**>
**>> dput(trans_matrix)
**>>
**>
**> structure(c(0.9426, 0.0047, 0, 4e-04, 0, 5e-04, 0, 0, 0, 0, 0,
**> 0, 9e-04, 0, 0, 0, 0, 0, 0.0308, 0.8205, 0.0254, 4e-04, 9e-04,
**> 0, 0, 0, 0, 4e-04, 0.0016, 0, 0, 0, 0, 0, 0, 0, 0.021, 0.1291,
**> 0.7978, 0.034, 0.0043, 0.0025, 0.0011, 3e-04, 9e-04, 4e-04, 0,
**> 0, 0, 0, 0, 0, 0, 0, 0.0056, 0.0394, 0.1174, 0.8366, 0.0509,
**> 0.0094, 0.0023, 0.0022, 0.0014, 0.0017, 0, 0, 0, 0, 0, 0, 0,
**> 0, 0, 0.0016, 0.0448, 0.0944, 0.8253, 0.0569, 0.011, 0.0051,
**> 6e-04, 0.0021, 0.0016, 0, 0, 0, 0, 0, 0, 0, 0, 0.0016, 0.0067,
**> 0.024, 0.0873, 0.8108, 0.0677, 0.0105, 0.004, 0.0017, 0, 9e-04,
**> 0, 0, 0.0039, 0.0045, 0, 0, 0, 0.0016, 0.0024, 0.0064, 0.0179,
**> 0.0833, 0.7838, 0.0758, 0.012, 0.0021, 0.0024, 0, 9e-04, 0, 0.0013,
**> 0.0015, 0, 0, 0, 0, 0.0024, 0.0014, 0.0068, 0.0202, 0.0908, 0.7694,
**> 0.0744, 0.0136, 0.0079, 0.0045, 9e-04, 0.0024, 0, 0, 0.0024,
**> 0, 0, 0.0016, 0, 0.0014, 0.004, 0.0089, 0.0175, 0.0983, 0.8066,
**> 0.1032, 0.0267, 0.0063, 0.0054, 0.0024, 0, 0.003, 0, 0, 0, 0,
**> 0.0024, 0, 0.0022, 0.002, 0.0079, 0.0178, 0.0632, 0.7605, 0.1422,
**> 0.0308, 0.0099, 0.006, 0.0013, 0.003, 0, 0, 0, 0, 0, 0, 0, 0.0017,
**> 0.0014, 0.0086, 0.0117, 0.0574, 0.6787, 0.1014, 0.0425, 0.006,
**> 0.0039, 0.003, 0.0047, 0, 0, 0, 0, 0, 3e-04, 0.0012, 0.0014,
**> 0.0022, 0.0109, 0.0227, 0.0589, 0.6814, 0.1058, 0.041, 0.0078,
**> 0.0045, 0.0024, 0, 0, 0, 6e-04, 0, 0, 2e-04, 0.0096, 0.0038,
**> 0.0034, 0.012, 0.0212, 0.0661, 0.6609, 0.159, 0.0362, 0.012,
**> 0.0024, 0, 0, 0, 0, 0, 0, 2e-04, 3e-04, 0.0016, 0.0029, 0.0058,
**> 0.0118, 0.029, 0.0597, 0.6205, 0.1332, 0.0331, 0.0047, 0, 0,
**> 0, 0, 0, 0, 0.001, 0.0011, 3e-04, 0.002, 0.0074, 0.0157, 0.0244,
**> 0.0515, 0.0916, 0.6546, 0.1203, 0.0353, 0, 0, 0, 0, 4e-04, 0,
**> 0, 3e-04, 0, 6e-04, 8e-04, 0.0047, 0.0199, 0.0172, 0.0325, 0.0699,
**> 0.7098, 0.1318, 0, 0, 0, 0, 0, 0, 5e-04, 0.0017, 0.0022, 0.004,
**> 0.0021, 0.0102, 0.0217, 0.0244, 0.0241, 0.044, 0.0737, 0.5271,
**> 0, 0, 0, 0, 7e-04, 0, 7e-04, 0.0023, 0.0019, 0.0014, 0.0062,
**> 0.0165, 0.0136, 0.0199, 0.0145, 0.044, 0.0316, 0.2894, 1), .Dim = c(18L,
**> 18L), .Dimnames = list(c("1", "2", "3", "4", "5", "6", "7", "8",
**> "9", "10", "11", "12", "13", "14", "15", "16", "17", "18"), c("AAA",
**> "AAp", "AA", "AAm", "Ap", "A", "Am", "BBBp", "BBB", "BBBm", "BBp",
**> "BB", "BBm", "Bp", "B", "Bm", "CCC.to.C", "D")))
**>
**>
**> ------ BEGIN PASTE ------
**>
**> # create a matrix of eigenvectors of the transition matrix
**> X <- eigen(trans_matrix)$vectors
**>
**> # create a diagonalized matrix of the eigenvalues of the transition matrix
**> L <- diag(eigen(trans_matrix)$values)
**>
**> # calculate inverse of matrix of eigenvectors of the transition matrix
**> X_inv <- solve(X)
**>
**> # calculate the 12th root of the eigenvalues in the diagonal matrix
**> L_star <- L ^ (1/12)
**>
**> # calculate the 12th root of the transition matrix
**> nth_root <- X_inv %*% L_star %*% X
**>
**> ------ END PASTE ------
**>
**>
**> References:
**>
**> [1] Hull, John. Risk Management and Financial Institutions. Prentice Hall, 2007.
**>
**> ______________________________________________
**> 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.
**>
*

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. Received on Fri 18 Jun 2010 - 10:21:47 GMT

Archive maintained by Robert King, hosted by
the discipline of
statistics at the
University of Newcastle,
Australia.

Archive generated by hypermail 2.2.0, at Fri 18 Jun 2010 - 11:20:32 GMT.

*
Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help.
Please read the posting
guide before posting to the list.
*