aboutsummaryrefslogtreecommitdiffstats
path: root/R/tranche_functions.R
diff options
context:
space:
mode:
Diffstat (limited to 'R/tranche_functions.R')
-rw-r--r--R/tranche_functions.R531
1 files changed, 531 insertions, 0 deletions
diff --git a/R/tranche_functions.R b/R/tranche_functions.R
new file mode 100644
index 00000000..0bccdf02
--- /dev/null
+++ b/R/tranche_functions.R
@@ -0,0 +1,531 @@
+library(lossdistrib)
+pos <- function(x){
+ pmax(x, 0)
+}
+
+trancheloss <- function(L, K1, K2){
+ pos(L - K1) - pos(L - K2)
+}
+
+trancherecov <- function(R, K1, K2){
+ pos(R - 1 + K2) - pos(R - 1 +K1)
+}
+
+tranche.cl <- function(L, R, cs, K1, K2, Ngrid=nrow(L), scaled=FALSE){
+ ## computes the couponleg of a tranche
+ ## if scaled is TRUE, scale it by the size of the tranche (K2-K1)
+ ## can make use of the fact that the loss and recov distribution are
+ ## truncated (in that case nrow(L) != Ngrid
+ if(K1==K2){
+ return( 0 )
+ }else{
+ support <- seq(0, 1, length=Ngrid)[1:nrow(L)]
+ size <- K2 - K1 - crossprod(trancheloss(support, K1, K2), L) -
+ crossprod(trancherecov(support, K1, K2), R)
+ sizeadj <- as.numeric(0.5 * (size + c(K2-K1, size[-length(size)])))
+ if(scaled){
+ return( 1/(K2-K1) * crossprod(sizeadj * cs$coupons, cs$df) )
+ }else{
+ return( crossprod(sizeadj * cs$coupons, cs$df) )
+ }
+ }
+}
+
+tranche.cl.scenarios <- function(l, r, cs, K1, K2, scaled=FALSE){
+ ## computes the couponleg of a tranche for one scenario
+ ## if scaled is TRUE, scale it by the size of the tranche (K2-K1)
+ ## can make use of the fact that the loss and recov distribution are
+ ## truncated (in that case nrow(L) != Ngrid
+ if(K1==K2){
+ return( 0 )
+ }else{
+ size <- K2 - K1 - trancheloss(l, K1, K2) - trancherecov(r, K1, K2)
+ sizeadj <- as.numeric(0.5 * (size + c(K2-K1, size[-length(size)])))
+ if(scaled){
+ return( 1/(K2-K1) * crossprod(sizeadj * cs$coupons, cs$df) )
+ }else{
+ return( crossprod(sizeadj * cs$coupons, cs$df) )
+ }
+ }
+}
+
+funded.tranche.pv <- function(L, R, cs, K1, K2, scaled = FALSE){
+ if(K1==K2){
+ return(0)
+ }else{
+ size <- K2 - K1 -trancheloss(L, K1, K2) - trancherecov(R, K1, K2)
+ sizeadj <- as.numeric(0.5 * (size + c(K2-K1, size[-length(size)])))
+ interest <- crossprod(sizeadj * cs$coupons, cs$df)
+ principal <- diff(c(0, trancherecov(R, K1, K2)))
+ principal[length(principal)] <- principal[length(principal)] + size[length(size)]
+ principal <- crossprod(cs$df, principal)
+ if(scaled){
+ pv <- (interest + principal)/(K2-K1)
+ }else{
+ pv <- (interest + principal)
+ }
+ return(pv)
+ }
+}
+
+tranche.pl <- function(L, cs, K1, K2, Ngrid=nrow(L), scaled=FALSE){
+ ## computes the protection leg of a tranche
+ ## if scaled
+ if(K1==K2){
+ return(0)
+ }else{
+ support <- seq(0, 1, length=Ngrid)[1:nrow(L)]
+ cf <- K2 - K1 - crossprod(trancheloss(support, K1, K2), L)
+ cf <- c(K2 - K1, cf)
+ if(scaled){
+ return( 1/(K2-K1) * crossprod(diff(cf), cs$df))
+ }else{
+ return( crossprod(diff(cf), cs$df))
+ }
+ }
+}
+
+tranche.pl.scenarios <- function(l, cs, K1, K2, scaled=FALSE){
+ ## computes the protection leg of a tranche
+ ## if scaled
+ if(K1==K2){
+ return(0)
+ }else{
+ cf <- K2 - K1 - trancheloss(l, K1, K2)
+ cf <- c(K2 - K1, cf)
+ if(scaled){
+ return( 1/(K2-K1) * as.numeric(crossprod(diff(cf), cs$df)))
+ }else{
+ return( as.numeric(crossprod(diff(cf), cs$df)))
+ }
+ }
+}
+
+tranche.pv <- function(L, R, cs, K1, K2, Ngrid=nrow(L)){
+ return( tranche.pl(L, cs, K1, K2, Ngrid) + tranche.cl(L, R, cs, K1, K2, Ngrid))
+}
+
+tranche.pv.scenarios <- function(l, r, cs, K1, K2){
+ return( tranche.pl.scenarios(l, cs, K1, K2, TRUE) +
+ tranche.cl.scenarios(l, r, cs, K1, K2, TRUE))
+}
+
+adjust.attachments <- function(K, losstodate, factor){
+ ## computes the attachments adjusted for losses
+ ## on current notional
+ return( pmin(pmax((K-losstodate)/factor, 0),1) )
+}
+
+BCtranche.legs <- function(index, K, rho, complement=FALSE){
+ ## computes the protection leg and couponleg of a 0-K tranche
+ ## if complement==TRUE, computes the protection leg and coupon leg of a K-1 tranche
+ if((K==0 && !complement) || (K==1 && complement)){
+ return(list(cl=0, pl=0))
+ }else if((K==1 && !complement) || (K==0 && complement)){
+ return(BCindex.pv(index))
+ }else{
+ dist <- BClossdistC(index$defaultprob, index$issuerweights, index$recov, rho, index$Z, index$w, index$N)
+ if(complement){
+ return(list(cl=tranche.cl(dist$L, dist$R, index$cs, K, 1),
+ pl=tranche.pl(dist$L, index$cs, K, 1)))
+ }else{
+ return(list(cl=tranche.cl(dist$L, dist$R, index$cs, 0, K),
+ pl=tranche.pl(dist$L, index$cs, 0, K)))
+ }
+ }
+}
+
+BCtranche.pv <- function(index, protection=FALSE, complement=FALSE){
+ ## computes the protection leg, couponleg, and bond price of a tranche
+ ## in the base correlation setting
+ ## if complement=FALSE compute the pvs starting from 0 (bottom-up skew)
+ ## if complement=TRUE compute the pvs starting from 1 (top-down skew)
+ pl <- rep(0, length(index$rho))
+ cl <- rep(0, length(index$rho))
+ for(i in seq_along(index$rho)){
+ temp <- BCtranche.legs(index, index$K[i], index$rho[i], complement)
+ pl[i] <- temp$pl
+ cl[i] <- temp$cl
+ }
+ dK <- diff(index$K)
+ plvec <- diff(pl)/dK
+ clvec <- diff(cl)/dK*index$tranches$running
+ if(complement){
+ plvec <- -plvec
+ clvec <- -clvec
+ }
+ if(protection){
+ bp <- -plvec-clvec
+ }else{
+ bp <- 1+plvec+clvec
+ }
+ return(list(pl=plvec, cl=clvec, bp=bp))
+}
+
+MFtranche.pv <- function(index, dist, protection=FALSE){
+ n.tranches <- length(index$K)-1
+ pl <- rep(0, n.tranches)
+ cl <- rep(0, n.tranches)
+ for(i in seq.int(n.tranches)){
+ pl[i] <- tranche.pl(dist$L, index$cs, index$K[i], index$K[i+1])
+ cl[i] <- tranche.cl(dist$L, dist$R, index$cs, index$K[i], index$K[i+1])
+ }
+ dK <- diff(index$K)
+ plvec <- pl/dK
+ clvec <- cl/dK*index$tranches$running
+ if(protection){
+ bp <- -plvec-clvec
+ }else{
+ bp <- 1+plvec+clvec
+ }
+ return(list(pl=plvec, cl=clvec, bp=bp))
+}
+
+adjust.skew <- function(index1, index2, method="ATM"){
+ #index1 is the index for which we already have computed the skew
+ #index2 is the index we're mapping to
+ # if method="ATM", do simple at the money mapping
+ # method="TLP", do tranche loss proportion mapping
+ # method="PM", do probability matching
+
+ K1 <- index1$K[-c(1,length(index1$K))]
+ K2 <- index2$K[-c(1,length(index2$K))]
+ aux <- function(x, index1, el1, skew, index2, el2, K2){
+ return(abs(ELtrunc(index1, x, skew(x))/el1-
+ ELtrunc(index2, K2, skew(x))/el2))
+ }
+ aux2 <- function(x, index1, skew, index2, K2){
+ return(abs(log(Probtrunc(index1, x, skew(x)))-
+ log(Probtrunc(index2, K2, skew(x)))))
+ }
+ if(method %in% c("ATM", "TLP")){
+ el1 <- EL(index1)
+ el2 <- EL(index2)
+ }
+ skew <- function(x){
+ #we cap the correlation at 0.99 and 0.01
+ f <- splinefun(K1, index1$rho[-c(1, length(index1$rho))], "natural")
+ return(pmax(pmin(f(x), 0.99), 0.01))
+ }
+ if(method=="ATM"){
+ K1eq <- el1/el2 * K2
+ }else if(method == "TLP"){
+ K1eq <- c()
+ m <- max(K2) + 0.3
+ for(K2val in K2){
+ prog <- optimize(aux, interval=c(0,m),
+ index1=index1, el1=el1, skew=skew,
+ index2=index2, el2=el2, K2=K2val)
+ K1eq <- c(K1eq, prog$minimum)
+ }
+ }else if (method=="PM"){
+ K1eq <- c()
+ m <- max(K2) + 0.25
+ for(K2val in K2){
+ prog <- optimize(aux2, interval=c(0, m),
+ index1=index1, skew=skew, index2=index2, K2=K2val)
+ K1eq <- c(K1eq, prog$minimum)
+ }
+ }
+ return(c(NA, skew(K1eq), NA))
+}
+
+theta.adjust.skew <- function(index, shortened=4, method="ATM"){
+ #ajust the correlation skew by doing ATM mapping on the expected loss
+ indexshort <- index
+ N <- nrow(index$cs)-shortened
+ indexshort$defaultprob <- indexshort$defaultprob[,1:N]
+ indexshort$cs <- indexshort$cs[1:N,]
+ return(adjust.skew(index, indexshort, method))
+}
+
+BCtranche.theta <- function(index, shortened=4, complement=FALSE, method="ATM"){
+ temp <- BCtranche.pv(index, complement=complement)
+ rho.adj <- theta.adjust.skew(index, shortened, method)
+ if(any(rho.adj[-c(1, length(rho.adj))]<=0)){
+ print("probable inverted skew: no adjustment")
+ }else{
+ index$rho <- rho.adj
+ }
+ N <- nrow(index$cs) - shortened
+ index$cs <- index$cs[1:N,]
+ index$defaultprob <- index$defaultprob[,1:N]
+ temp2 <- BCtranche.pv(index, complement=complement)
+ temp3 <- BCtranche.delta(index, complement=complement)
+ return(data.frame(theta=temp2$bp-temp$bp+index$tranches$running,
+ fw.delta=temp3$delta))
+}
+
+BCtranche.delta <- function(index, complement=FALSE){
+ ## computes the tranche delta (on current notional) by doing a proportional
+ ## blip of all the curves
+ ## if complement is False, then computes deltas bottom-up
+ ## if complement is True, then computes deltas top-down
+ eps <- 1e-4
+ index$N <- 301 ## for gamma computations we need all the precision we can get
+ ## we build a lit of 4 indices with various shocks
+ index.list <- lapply(c(0, eps, -eps, 2*eps), function(x){
+ if(x==0){
+ return(index)
+ }else{
+ newindex <- index
+ newindex$portfolio <- tweakportfolio(newindex$portfolio, x)
+ newindex$defaultprob <- 1 - SPmatrix(newindex$portfolio, length(index$cs$dates))
+ return(newindex)
+ }
+ })
+ bp <- matrix(0, length(index$K)-1, length(index.list))
+ indexbp <- rep(0, length(index.list))
+ for(j in seq_along(index.list)){
+ indexbp[j] <- BCindex.pv(index.list[[j]])$bp
+ bp[,j] <- BCtranche.pv(index.list[[j]], complement=complement)$bp
+ }
+
+ deltas <- (bp[,2]-bp[,3])/(indexbp[2]-indexbp[3])*tranche.factor(index)/index$factor
+ deltasplus <- (bp[,4]-bp[,1])/(indexbp[4]-indexbp[1])*tranche.factor(index)/index$factor
+ gammas <- (deltasplus-deltas)/(indexbp[2]-indexbp[1])/100
+
+ return( data.frame(delta=deltas, gamma=gammas) )
+}
+
+MFtranche.delta <- function(index){
+ ## computes the tranche delta (on current notional) by doing a proportional
+ ## blip of all the curves
+ ## if complement is False, then computes deltas bottom-up
+ ## if complement is True, then computes deltas top-down
+ eps <- 1e-4
+ index$Ngrid <- 301 ## for gamma computations we need all the precision we can get
+ ## we build a lit of 4 indices with various shocks
+ index.list <- lapply(c(0, eps, -eps, 2*eps), function(x){
+ if(x==0){
+ return(index)
+ }else{
+ newindex <- index
+ newindex$portfolio <- tweakportfolio(newindex$portfolio, x)
+ newindex$defaultprob <- 1 - SPmatrix(newindex$portfolio, length(index$cs$dates))
+ return(newindex)
+ }
+ })
+ bp <- matrix(0, length(index$K)-1, length(index.list))
+ indexbp <- rep(0, length(index.list))
+ for(j in seq_along(index.list)){
+ indexbp[j] <- BCindex.pv(index.list[[j]])$bp
+ dist <- MFdist(index.list[[j]])
+ bp[,j] <- BCtranche.pv(index.list[[j]], dist)
+ }
+
+ deltas <- (bp[,2]-bp[,3])/(indexbp[2]-indexbp[3])*tranche.factor(index)/index$factor
+ deltasplus <- (bp[,4]-bp[,1])/(indexbp[4]-indexbp[1])*tranche.factor(index)/index$factor
+ gammas <- (deltasplus-deltas)/(indexbp[2]-indexbp[1])/100
+
+ return( list(deltas=deltas, gammas=gammas) )
+}
+
+BCtranche.corr01 <- function(index, eps=0.01, complement=FALSE){
+ ##does a parallel shift of the skew and computes the change in pv
+ before <- BCtranche.pv(index, complement=complement)
+ index$rho[-1] <- index$rho[-1]+eps
+ after <- BCtranche.pv(index, complement=complement)
+ return(after$bp-before$bp)
+}
+
+EL <- function(index, discounted=TRUE, shortened=0){
+ ## computes the expected loss of a portfolio (time discounted if discounted is TRUE)
+ ## given the default curves and recovery
+ ## should be very close to the protection leg of the portfolio of cds
+ ## index should be a list with issuerweights, recov, defaultprob and cs parameters
+ ## shortened: number of quarters to shorten the maturity by
+ Ncol <- ncol(index$defaultprob)-shortened
+ ELvec <- as.numeric(crossprod(index$issuerweights * (1-index$recov), index$defaultprob[,1:Ncol]))
+ if(!discounted){
+ return( ELvec[length(ELvec)] )
+ }else{
+ return( sum(index$cs$df[1:Ncol]*diff(c(0, ELvec))) )
+ }
+}
+
+BCindex.pv <- function(index, discounted=TRUE, shortened=0){
+ Ncol <- ncol(index$defaultprob)-shortened
+ ELvec <- as.numeric(crossprod(index$issuerweights * (1-index$recov), index$defaultprob[,1:Ncol]))
+ size <- 1-as.numeric(crossprod(index$issuerweights, index$defaultprob[,1:Ncol]))
+ sizeadj <- 0.5*(c(1, size[-length(size)])+size)
+ if(!discounted){
+ pl <- -ELvec[length(ELvec)]
+ cl <- as.numeric(crossprod(index$cs$coupons[1:Ncol], sizeadj))
+ bp <- 1+cl+pl
+ }else{
+ pl <- -sum(index$cs$df[1:Ncol]* diff(c(0, ELvec)))
+ cl <- as.numeric(crossprod(index$cs$coupons[1:Ncol], sizeadj * index$cs$df[1:Ncol] ))
+ bp <- 1+cl+pl
+ }
+ bp <- 1+cl*index$quotes$spread+pl
+ return(list(pl=pl, cl=cl, bp=bp))
+}
+
+ELtrunc <- function(index, K, rho){
+ ## computes the expected loss of a portfolio below strike K
+ ## could be written faster by using a truncated version of lossdist
+ ## index should be a list with issuerweights, recov, defaultprob and cs parameters
+ dist <- BClossdistC(index$defaultprob, index$issuerweights, index$recov, rho, index$Z, index$w, index$N)
+ return( -tranche.pl(dist$L, index$cs, 0, K))
+}
+
+Probtrunc <- function(index, K, rho){
+ dist <- BClossdistC(index$defaultprob, index$issuerweights, index$recov, rho, index$Z, index$w, index$N)
+ p <- cumsum(dist$L[,ncol(dist$L)])
+ support <- seq(0, 1, length=index$N)
+ probfun <- splinefun(support, p, method="hyman")
+ return(probfun(K))
+}
+
+
+BCstrikes <- function(index, K, rho) {
+ ## computes the strikes as a percentage of expected loss
+ ## Kmodified is the current attachment points (adjusted for losses)
+ el <- EL(index)
+ ELvec <- c()
+ for(i in 2:length(K)){
+ ELvec <- c(ELvec, ELtrunc(index, K[i], rho[i]))
+ }
+ return( ELvec/el )
+}
+
+tranche.factor <- function(index){
+ ## compute the factor to convert from delta on current notional to delta on original notional
+ ## K1 and K2 original strikes
+ return( diff(index$K)/diff(index$K.orig)*index$factor )
+}
+
+MFupdate.prob <- function(Z, w, rho, defaultprob, useC = TRUE){
+ ## update the probabilites based on a non gaussian factor
+ ## distribution so that the pv of the cds stays the same.
+ p <- matrix(0, nrow(defaultprob), ncol(defaultprob))
+ fit.prob <- if(useC) fit.probC else fit.prob
+ for(i in 1:nrow(defaultprob)){
+ for(j in 1:ncol(defaultprob)){
+ p[i,j] <- fit.prob(Z, w, rho[i], defaultprob[i,j])
+ }
+ }
+ return( p )
+}
+
+MFlossrecovdist.prepay <- function(w, Z, rho, defaultprob, defaultprobmod, prepayprob, prepayprobmod,
+ issuerweights, recov, Ngrid=2*length(issuerweights)+1, defaultflag=FALSE){
+ ## computes the loss and recovery distribution using the modified factor distribution
+ n.credit <- length(issuerweights)
+ n.int <- length(w)
+ Rstoch <- array(0, dim=c(n.int, n.credit, ncol(defaultprob)))
+ for(t in 1:ncol(defaultprob)){
+ for(i in 1:n.credit){
+ Rstoch[,i,t] <- stochasticrecov(recov[i], 0, Z, w, rho, defaultprob[i,t], defaultprobmod[i,t])
+ }
+ }
+ parf <- function(i){
+ dpshocked <- apply(defaultprobmod, 2, shockprob, rho=rho, Z=Z[i])
+ ppshocked <- apply(prepayprobmod, 2, shockprob, rho=rho, Z=-Z[i])
+ S <- 1 - Rstoch[i,,]
+ dist <- lossrecovdist.term(dpshocked, ppshocked, issuerweights, S, Ngrid, defaultflag)
+ }
+ L <- matrix(0, Ngrid, ncol(defaultprob))
+ R <- matrix(0, Ngrid, ncol(defaultprob))
+ for(i in 1:length(w)){
+ dist <- parf(i)
+ L <- L + dist$L * w[i]
+ R <- R + dist$R * w[i]
+ }
+ return( list(L=L, R=R) )
+}
+
+MFlossdist.joint <- function(cl, w, Z, rho, defaultprob, defaultprobmod, issuerweights, recov,
+ Ngrid=2*length(issuerweights)+1, defaultflag=FALSE){
+ ## rowSums(Q) is the default/loss distribution depending if
+ ## defaultflag is TRUE or FALSE (default setting is FALSE)
+ ## colSums(Q) is the recovery distribution
+ ## so that recovery is the y axis and L/D is the x axis
+ ## if we use the persp function, losses is the axis facing us,
+ ## and R is the axis going away from us.
+ n.credit <- length(issuerweights)
+ n.int <- lenth(w)
+ Rstoch <- array(0, dim=c(n.int, n.credit, ncol(defaultprob)))
+ for(t in 1:ncol(defaultprob)){
+ for(i in 1:n.credit){
+ Rstoch[,i,t] <- stochasticrecov(recov[i], 0, Z, w, rho, defaultprob[i,t], defaultprobmod[i,t])
+ }
+ }
+ parf <- function(i){
+ pshocked <- apply(defaultprobmod, 2, shockprob, rho=rho, Z=Z[i])
+ S <- 1 - Rstoch[i,,]
+ dist <- lossrecovdist.joint.term(pshocked, 0, issuerweights, S, Ngrid, defaultflag)
+ gc()
+ return(dist)
+ }
+ temp <- parSapply(cl, 1:length(w), parf)
+ clusterCall(cl, gc)
+ Q <- array(0, dim=c(ncol(defaultprob), Ngrid, Ngrid))
+ for(i in 1:length(w)){
+ Q <- Q + w[i]*array(temp[,i], dim=c(ncol(defaultprob), Ngrid, Ngrid))
+ }
+ return( Q )
+}
+
+MFlossdist.prepay.joint <- function(w, Z, rho, defaultprob, defaultprobmod,
+ prepayprob, prepayprobmod, issuerweights, recov,
+ Ngrid=2*length(issuerweights)+1, defaultflag=FALSE){
+ ## rowSums is the loss distribution
+ ## colSums is the recovery distribution
+ ## so that recovery is the y axis and L is the x axis
+ ## if we use the persp function, losses is the axis facing us,
+ ## and R is the axis going away from us.
+ n.credit <- length(issuerweights)
+ n.int <- length(w)
+ Rstoch <- array(0, dim=c(n.credit, n.int, ncol(defaultprob)))
+
+ for(t in 1:ncol(defaultprob)){
+ for(i in 1:n.credit){
+ Rstoch[i,,t] <- stochasticrecovC(recov[i], 0, Z, w, rho[i],
+ defaultprob[i,t], defaultprobmod[i,t])
+ }
+ }
+ Q <- array(0, dim=c(ncol(defaultprob), Ngrid, Ngrid))
+ for(t in 1:ncol(defaultprob)){
+ S <- 1 - Rstoch[,,t]
+ Q[t,,] <- lossdistC.prepay.jointZ(defaultprobmod[,t], prepayprobmod[,t],
+ issuerweights, S, Ngrid, defaultflag, rho, Z, w)
+ }
+ return( Q )
+}
+
+MFrecovery <- function(index, defaultprobmod){
+ n.credit <- length(index$issuerweights)
+ n.int <- length(index$Z)
+ Rstoch <- array(0, dim=c(n.credit, n.int, ncol(index$defaultprob)))
+ rho <- rep(0.45, n.credit)
+ for(t in 1:ncol(index$defaultprob)){
+ for(i in 1:n.credit){
+ Rstoch[i,,t] <- stochasticrecovC(index$recov[i], 0, index$Z, index$w.mod,
+ rho[i], index$defaultprob[i,t], defaultprobmod[i,t])
+ }
+ }
+ return( Rstoch )
+}
+
+MFlossdist <- function(index){
+ n.credit <- length(index$issuerweights)
+ rho <- rep(0.45, n.credit)
+ defaultprobmod <- MFupdate.probC(index$Z, index$w.mod, rho, index$defaultprob)
+ n.credit <- length(index$issuerweights)
+ n.int <- length(index$Z)
+ Rstoch <- MFrecovery(index, defaultprobmod)
+ Lw <- matrix(0, index$N, n.int)
+ Rw <- matrix(0, index$N, n.int)
+ L <- matrix(0, index$N, ncol(index$defaultprob))
+ R <- matrix(0, index$N, ncol(index$defaultprob))
+ for(t in 1:ncol(index$defaultprob)){
+ S <- 1 - Rstoch[,,t]
+ Lw <- lossdistCZ(defaultprobmod[,t], index$issuerweights, S, index$N, 0, rho, index$Z)
+ Rw <- lossdistCZ(defaultprobmod[,t], index$issuerweights, 1-S, Ngrid, 0, rho, index$Z)
+ L[,t] <- Lw%*%index$w.mod
+ R[,t] <- Rw%*%index$w.mod
+ }
+ return(list(L=L, R=R))
+}
generated by cgit v1.2.3-70-g09d2 (git 2.52.0) at 2026-01-10 17:54:57 +0000