How can I compute the covariance of bond returns for a universe of bonds? Unlike equities, bonds have a finite maturity and thus there is a price path dependence as the bond matures and expires at par value. In the case of equities, one can use historical returns to calculate an empirical covariance matrix, a statistical risk modeling approach using PCA, or use some fundamental factor model, where I would estimate "factor returns" based on a linear regression of past equity returns against their exposure to said factors. With bonds, all of these methods do not easily translate because the returns used are non-stationary as their volatility shrinks as maturity approaches. This is where I get stuck.
I have seen people create structural models using a Vasicek model (or a Cox-Ingersoll-Ross...CIR) model with a default probability included to get a volatility estimate of a bond trading. Meaning, if I have a 6-month bond and it was issued today I can use these models to derive an estimate of its return volatility today without having actually seen any real returns yet. However, how can I arrive at a covariance between this hypothetical 6-month bond and another 5-year bond from a different company?
I have also seen simulation approaches, where people have suggested re-sampling the historical forward curve to create synthetic yield curves and then calculate the "return" of the bond as its new price, after re-pricing the bond with the new simulated yield curve, divided by the actual market price. This does allow me to obtain a vector of returns that I can calculate a covariance matrix for, but it seems computationally intensive.
Any thoughts or pointers on things to read on this topic would be appreciated. Thanks!