upper_triangular_to_ascending_triangle = function(M){
    ##
    ## Converts a matrix from the upper triangular form used for most of the book
    ## into the ascending triangle form used in the chapter "Interest Rate Derivatives".
    ##
    ## In the first form (the upper triangular form) an
    ## "up" step is on the same horizontal row and a 
    ## "down" step is one row down from the current row.
    ##
    ## In the second form (the ascending triangle form) an
    ## "up" step is on one row above the curent row and a 
    ## "down" step is on the same row as the current row.
    ##
    m = dim(M)[1]
    n = dim(M)[2]
    M_out = matrix(data=rep(NA, m*n), nrow=m, ncol=n, byrow=TRUE)
    for( jj in 1:n ){
        for( ii_target in 1:jj ){
            ii_destination = m-jj+ii_target
            M_out[ii_destination, jj] = M[ii_target, jj]
        }
    }
    rownames(M_out) = rownames(M)
    colnames(M_out) = colnames(M)
    return(M_out)
}


short_rate_lattice_to_elementary_prices = function(rrate){
    ##
    ## okular ../EBook/Investment_Science.djvu -p 412 & # the forward equations 
    ##
    n_rows = dim(rrate)[1]
    n_cols = dim(rrate)[2] + 1 
    n_years = n_cols - 1

    P0 = matrix(NA, nrow=n_rows, ncol=n_cols)
    P0[n_rows, 1] = 1.0
    for( kk in seq(0, n_cols-2) ){ ## the physical time index
        for( ss in 0:(kk+1) ){ ## physical up/down index

            if( ss==0 ){
                ii = n_rows ## the array indices where we came from 
                jj = kk + 1
                d_k_0 = 1/(1+rrate[ii, jj])
                P0[ii, jj+1] = 0.5 * d_k_0 * P0[ii, jj]
            }else if( ss==(kk+1) ){
                ii = n_rows - ss + 1 ## the array indices where we came from 
                jj = kk+1
                d_k_k = 1/(1+rrate[ii, jj])
                P0[ii-1, jj+1] = 0.5 * d_k_k * P0[ii, jj]
            }else{ 
                ii = n_rows - ss + 1 
                jj = kk+1
                d_k_sm1 = 1/(1 + rrate[ii, jj])
                P0[ii-1, jj+1] = 0.5 * d_k_sm1 * P0[ii, jj]

                ii = n_rows - ss
                jj = kk+1
                d_k_s = 1/(1 + rrate[ii, jj])
                P0[ii, jj+1] = P0[ii, jj+1] + 0.5 * d_k_s * P0[ii, jj]
            }

        }
    }
    colnames(P0) = paste0('t_', seq(0, n_years))
    P0
}


constant_b_ho_lee_model_to_short_rate = function(ak, b, n_years){
    ##
    ## From the Ho-Lee model: 
    ##
    ## r_{ks} = a_k + b_k s 
    ##
    ## for a constant b_k i.e. b_k = b
    ##
    ## we compute the short rates lattice
    ##
    n_rows = n_years + 1
    n_cols = n_years + 1 

    rrate = matrix(data=NA, nrow=n_rows, ncol=n_cols)
    for( kk in seq(0, n_years) ){ ## the physical time index 
        for( ss in 0:kk ){ ## physical up/down index
            ii = n_rows - ss
            jj = kk+1
            rrate[ii, jj] = ak[jj] + b * ss
	}
    }
    rrate * 1.e-2 ## convert into numbers (not percents)
}


constant_b_black_derman_toy_model_to_short_rate = function(ak, b, n_years){
    ##
    ## From the Black-Derman-Toy model: 
    ##
    ## r_{ks} = a_k e^{b_k s} 
    ##
    ## for a constant b_k i.e. b_k = b
    ##
    ## we compute the short rates lattice
    ##
    n_rows = n_years + 1
    n_cols = n_years + 1 

    rrate = matrix(data=NA, nrow=n_rows, ncol=n_cols)
    for( kk in seq(0, n_years) ){ ## the physical time index
        for( ss in 0:kk ){ ## physical up/down index
            ii = n_rows - ss
            jj = kk+1
            rrate[ii, jj] = ak[jj] * exp(b * ss)
	}
    }
    rrate * 1.e-2 ## convert into numbers (not percents)
}


model_term_structure_match_min_fn = function(spot_rates, ak, b, model_to_use='ho_lee'){ 
    ##
    ##
    ##
    n_years = length(spot_rates)-1

    if( model_to_use=='ho_lee' )
        rrate = constant_b_ho_lee_model_to_short_rate(ak, b, n_years)
    else if( model_to_use=='black_derman_toy' ){
        rrate = constant_b_black_derman_toy_model_to_short_rate(ak, b, n_years)
    }

    E0 = short_rate_lattice_to_elementary_prices(rrate)

    E0[is.na(E0)] = 0 ## sum the columns of the elemntary price grid
    E00 = colSums(E0)

    ks = 0:(n_years+1)
    sk = (1/E00)^(1/ks) - 1
    sk_pct = sk[-1] * 100 ## drop the first element and conver to a percent

    ##mean(abs(spot_rates - sk_pct), na.rm=TRUE)
    sum((spot_rates - sk_pct)^2)
}