# This file contains functions to estimate the goodness of fit of the death models
##########################################################################################
# einfaches Modell
GoodnessEinfachData = function(){
# Beschreibung: Diese Funktion bestimmt den Fehler, der auf dem original Datensatz
# mit der einfachen Schätzmethode erzeugt wurde.
# select half the data points
firsthalf = min(deathrates1879westmatrix[,1]) +
length(min(deathrates1879westmatrix[,1]):max(deathrates1879westmatrix[,1]))/2
rates = subset(deathrates1879westmatrix, deathrates1879westmatrix[,1] <= firsthalf)
exposure = subset(exposure1879westmatrix, exposure1879westmatrix[,1] <= firsthalf)
erg = einfachstesModellAlterProjektion(rates, exposure)
X = erg$X
Y_first_half = erg$Y
# estimate the parameters
par = einfachstesModellAlterEstimateParameters(X, Y_first_half)
# select the other half of the data points
rates_last = subset(deathrates1879westmatrix, deathrates1879westmatrix[,1] > firsthalf)
exposure_last = subset(exposure1879westmatrix, exposure1879westmatrix[,1] > firsthalf)
erg_last = einfachstesModellAlterProjektion(rates_last, exposure_last)
Y_second_half = erg_last$Y
# just to check
pdf("../../1 Doku/graphics/EinfachData.pdf", width = 10, height = 8)
plot(X, (2*pi*par$sigma^2)^(-1/2) * exp(- (X-par$mu)^2/(2*par$sigma^2)), type = "l",
ylim = c(0,0.04), ylab = "Sterbewahrscheinlichkeit", xlab = "Alter")
points(X, Y_first_half, pch = 1)
points(X,Y_second_half, pch = 4)
legend( x="topleft", cex = 1.5,
legend=c("Datensatz bis 1932","Datensatz ab 1932"), lwd=1, lty=c(0,0), pch=c(1,4))
dev.off()
# calculate the errors
errors = (Y_second_half - (2*pi*par$sigma^2)^(-1/2) * exp(- (X-par$mu)^2/(2*par$sigma^2)))^2
# estimate the error of the whole data set
erg_full = einfachstesModellAlterProjektion(deathrates1879westmatrix,
exposure1879westmatrix)
X_full = erg_full$X
Y_full = erg_full$Y
par_full = einfachstesModellAlterEstimateParameters(X_full, Y_full)
# estimate the variance
sigma_epsilon_hat = sqrt(sum((Y_full - (2*pi*par_full$sigma^2)^(-1/2) *
exp(- (X_full-par_full$mu)^2/(2*par_full$sigma^2)))^2)/
(length(Y_full)))
# estimate the error
return(sum(errors))
}
# estimate parameters einfaches modell
GoodnessEinfachSim = function(){
# Beschreibung: Diese Funktion bestimmt den Fehler, der auf dem simulierten Datensatz
# mit der einfachen Schätzmethode erzeugt wurde.
# start by generating a plot of the historical data and the generated data
pdf("../../1 Doku/graphics/SampleDataSimple.pdf", width = 10, height = 8)
erg = einfachstesModellAlterProjektion(deathrates1879westmatrix, exposure1879westmatrix)
X = erg$X
Y = erg$Y
erg = einfachstesModellAlterEstimateParameters(X, Y)
plot(simple_data[,1], simple_data[,3], pch = 4, xlab = "Alter", ylab = "Sterblichkeit")
lines(simple_data[,1], simple_data[,2])
curve((2*pi*erg$sigma^2)^(-1/2) * exp(- (x-erg$mu)^2/(2*erg$sigma^2)), add = T,
lty = "dashed")
legend("topleft", legend=c("errorless", "original"),
col=c("black", "black"), lty=1:2, cex=1.5)
dev.off()
# initialize constants
m = 10000
X = 0:95
# storage space for parameter estimators and diviation from true parameter
ms = rep(NA, m)
ss = rep(NA, m)
merror = rep(NA, m)
serror = rep(NA, m)
errors = rep(NA, m)
# estimate the errors
for(i in (1+2):(m+2)){
par = einfachstesModellAlterEstimateParameters(simple_data[,1], simple_data[,i])
ms[i-2] = par$mu
ss[i-2] = par$sigma
merror[i-2] = (par$mu- mu_data)^2
serror[i-2] = (par$sigma - sigma_data)^2
errors[i-2] = (par$mu- mu_data)^2 + (par$sigma - sigma_data)^2
}
# plot the different errors
pdf("../../1 Doku/graphics/ErrorSimpleDetailed.pdf", width = 10, height = 8)
par(mfrow=c(1,2))
hist(ms)
hist(ss)
dev.off()
par(mfrow=c(1,1))
return(mean(errors))
}
##########################################################################################
# The Lee-Carter Modell
GoodnessLeeData = function(){
# using the first half of the dataset estimate the second half and calculate the errors
# estimate the parameters
firsthalf = min(deathrates1965west[,1]) +
length(min(deathrates1965west[,1]):max(deathrates1965west[,1]))/2
Geburtsjahre = min(deathrates1965west[,1]):firsthalf
Alter = min(deathrates1965west[,2]):95
Deathrates = subset(deathrates1965west, deathrates1965west[,1] <= firsthalf)
erg = Lee(min(deathrates1965west[,1]):firsthalf, Alter,
subset(Deathrates, Deathrates[,1] <= firsthalf))
alpha_a = erg$alpha_a
beta_a = erg$beta_a
gamma_t = erg$gamma_t
nu = erg$nu
# calculate sigma xi for the simulations
sigma_xi_vec = rep(NA, length(Geburtsjahre)-1)
for(i in 1:(length(Geburtsjahre)-1)){
sigma_xi_vec[i] = (gamma_t[i+1] - gamma_t[i] - nu)^2
}
sigma_xi = sqrt((1/(length(Geburtsjahre)))*sum(sigma_xi_vec))
# predict
YearsToPredict = max(deathrates1965west[,1]) - firsthalf
nu_vec = seq(from = nu, to = nu * YearsToPredict, by = nu)
gamma_vec = c(gamma_t, gamma_t[length(gamma_t)] + nu_vec)
# generate the alpha matrix and beta matrix -> generate EstRates matrix
ones = matrix(rep(1, 1 + max(deathrates1965west[,1]) - min(deathrates1965west[,1])),
nrow = 1)
alpha_matrix = matrix(alpha_a, ncol = 1) %*% ones
beta_vec = matrix(beta_a, nrow = 1)
beta_gamma = t(beta_vec) %*% t(gamma_vec)
EstRates = exp(alpha_matrix + beta_gamma)
#EstRates = alpha_matrix + beta_gamma #for error estimation
# save the data for the simulation
alpha_data = alpha_a
beta_data = beta_a
nu_data = nu
gamma_historical = gamma_t
devtools::use_data(alpha_data, overwrite = T)
devtools::use_data(beta_data, overwrite = T)
devtools::use_data(nu_data, overwrite = T)
devtools::use_data(gamma_historical, overwrite = T)
plot(0:31, gamma_historical)
#lines(4:35, gamma_historical, lty = "dotted")
# plot the gamma parameter of the historical data set
pdf("../../1 Doku/graphics/gamma_historical.pdf", width = 10, height = 8)
plot(0:31, gamma_historical, type = "l", ylab = "alpha")
lines(0:31,0:31*nu+max(gamma_historical), lty = "dashed")
lines(0:31, rep(0,32), lty = "dotdash")
dev.off()
# calculate the errors
errors = rep(NA, ((1+max(deathrates1965west[,1])) - min(deathrates1965west[,1]))/2)
Geburtsjahre_all=max(deathrates1965west[,1]) : firsthalf # The last year comes first
j=1
for(i in Geburtsjahre_all){
errors[j] = sum((EstRates[,i - 1955] -
subset(deathrates1965west[,3], deathrates1965west[,1] == i))^2)
#log(subset(Deathrates[,3], Deathrates[,1] == i)))^2)#for error estimation
j=j+1
}
sd_rates = sum(errors)/length(errors)
# plot different examples of death probability estimation
pdf("../../1 Doku/graphics/ErrorLeeData.pdf", width = 10, height = 8)
par(mfrow = c(2,2))
# first year:1956
plot(Alter, subset(deathrates1965west[,3], deathrates1965west[,1] == 1956), type = "l",
ylab = "Sterblichkeit", main = "1956")
lines(Alter, EstRates[,1956-1955], lty = "dashed")
# end of first half:1987
plot(Alter, subset(deathrates1965west[,3], deathrates1965west[,1] == 1987), type = "l",
ylab = "Sterblichkeit", main = "1987")
lines(Alter, EstRates[,1987-1955], lty = "dashed")
# middle of prediction:2002
plot(Alter, subset(deathrates1965west[,3], deathrates1965west[,1] == 2002), type = "l",
ylab = "Sterblichkeit", main = "2002")
lines(Alter, EstRates[,2002-1955], lty = "dashed")
# end of prediction:2017
plot(Alter, subset(deathrates1965west[,3], deathrates1965west[,1] == 2017), type = "l",
ylab = "Sterblichkeit", main = "2017")
lines(Alter, EstRates[,2017-1955], lty = "dashed")
par(mfrow = c(1,1))
dev.off()
# plot different examples of death probability estimation
# log plot
pdf("../../1 Doku/graphics/LogErrorLeeData.pdf", width = 10, height = 8)
par(mfrow = c(2,2))
# first year:1956
plot(Alter, log(subset(deathrates1965west[,3], deathrates1965west[,1] == 1956)), type = "l",
ylab = "Sterblichkeit", main = "1956")
lines(Alter, log(EstRates[,1956-1955]), lty = "dashed")
# end of first half:1987
plot(Alter, log(subset(deathrates1965west[,3], deathrates1965west[,1] == 1987)), type = "l",
ylab = "Sterblichkeit", main = "1987")
lines(Alter, log(EstRates[,1987-1955]), lty = "dashed")
# middle of prediction:2002
plot(Alter, log(subset(deathrates1965west[,3], deathrates1965west[,1] == 2002)), type = "l",
ylab = "Sterblichkeit", main = "2002")
lines(Alter, log(EstRates[,2002-1955]), lty = "dashed")
# end of prediction:2017
plot(Alter, log(subset(deathrates1965west[,3], deathrates1965west[,1] == 2017)), type = "l",
ylab = "Sterblichkeit", main = "2017")
lines(Alter, log(EstRates[,2017-1955]), lty = "dashed")
par(mfrow = c(1,1))
dev.off()
return(mean(errors))
}
GoodneesLeeSim = function(){
# estimate parameters and errors on the whole data set
Zeitraum = min(complex_period_data[,1]):max(complex_period_data[,1])
Alter = min(complex_period_data[,2]):max(complex_period_data[,2])
est_data = rep(NA, length(Alter)*length(Zeitraum))
m = 10000
alphaerror = rep(NA, m)
betaerror = rep(NA, m)
gammaerror = rep(NA, m)
nuerror = rep(NA, m)
errors = rep(NA, m)
for(i in 1:m){
erg = Lee(Zeitraum, Alter, cbind(complex_period_data[,1:2], complex_period_data[,3+i]))
alphaerror[i] = sum((erg$alpha - alpha_data)^2)
betaerror[i] = sum((erg$beta - beta_data)^2)
gammaerror[i] = sum((erg$gamma - gamma_data)^2)
nuerror[i] = sum((erg$nu - nu_data)^2)
errors[i] = sum((erg$alpha - alpha_data)^2) + sum((erg$beta - beta_data)^2) +
sum((erg$gamma - gamma_data)^2) + sum((erg$nu - nu_data)^2)
}
pdf("../../1 Doku/graphics/ErrorLeeDetailed.pdf", width = 10, height = 8)
par(mfrow=c(2,2))
hist(alphaerror)
hist(betaerror)
hist(gammaerror)
hist(nuerror)
par(mfrow=c(1,1))
dev.off()
error_per_parameter = mean(errors)/(2*95+40+1)
return(errors = mean(errors))
}
GoodnessLeeSimPredict = function(){
# estimate the goodness of the Lee-Carter Model as predictor; to this end use the first
# half of the simulated data set to predict the second half and compare to the true values
firsthalf = round(min(complex_period_data[,1]) +
length(min(complex_period_data[,1]):max(complex_period_data[,1]))/2) # Zeitraum 1,...,40
Zeitraum = min(complex_period_data[,1]):firsthalf
Zeitraum_full = min(complex_period_data[,1]):max(complex_period_data[,1])
Alter = min(complex_period_data[,2]):max(complex_period_data[,2])
# calculate the mean of EstRates
m = 100
EstRatessum = matrix(rep(0, length(Alter)*length(Zeitraum_full)), ncol = 41)
for(i in 1:m){
Deathrates = cbind(subset(complex_period_data[,1:2], complex_period_data[,1] <= firsthalf),
subset(complex_period_data[,3+i], complex_period_data[,1] <= firsthalf))
erg = Lee(Zeitraum, Alter, Deathrates)
ones = matrix(rep(1, 1 + max(complex_period_data[,1]) - min(complex_period_data[,1])),
nrow = 1)
alpha_matrix = matrix(erg$alpha_a, ncol = 1) %*% ones
beta_vec = matrix(erg$beta_a, nrow = 1)
YearsToPredict = max(complex_period_data[,1]) - firsthalf
gamma_t = erg$gamma_t
nu = erg$nu
nu_vec = seq(from = nu, to = nu * YearsToPredict, by = nu)
gamma_vec = c(gamma_t, gamma_t[length(gamma_t)] + nu_vec)
beta_gamma = t(beta_vec) %*% t(gamma_vec)
EstRates = exp(alpha_matrix + beta_gamma)
EstRatessum = EstRatessum + EstRates
}
# now take the mean
EstRates = EstRatessum/m
# calculate errors
errors = rep(NA, (1+max(complex_period_data[,1])))
for(k in 1:41){
errors[k] = sum((EstRates[,k] - complex_period_data[(1:96)+96*(k-1),3])^2)
}
# estimate the death probabilities on the whole data set
EstRates_full_sum = matrix(rep(0, length(Alter)*length(Zeitraum_full)), ncol = 41)
for(i in 1:m){
est_data_full = Lee(Zeitraum_full, Alter, cbind(complex_period_data[,1:2],
complex_period_data[,3+i]))
ones_full = matrix(rep(1, 1 + max(complex_period_data[,1]) - min(complex_period_data[,1])),
nrow = 1)
alpha_matrix_full = matrix(est_data_full$alpha_a, ncol = 1) %*% ones
beta_vec_full = matrix(est_data_full$beta_a, nrow = 1)
gamma_vec_full= matrix(est_data_full$gamma_t, ncol = 1)
beta_gamma_full = t(beta_vec_full) %*% t(gamma_vec_full)
EstRates_full = exp(alpha_matrix_full + beta_gamma_full)
EstRates_full_sum = EstRates_full_sum + EstRates_full
}
# now take the mean
EstRates_full = EstRates_full_sum/m
# calculate errors
errors_full = rep(NA, (1+max(complex_period_data[,1])))
for(k in 1:41){
errors_full[k] = sum((EstRates_full[,k] - complex_period_data[(1:96)+96*(k-1),3])^2)
}
# plot the parameters
pdf("../../1 Doku/graphics/ParameterLee.pdf", width = 10, height = 8)
par(mfrow=c(2,2))
plot(Alter, alpha_data, type = "l", ylab = "alpha")
lines(Alter, est_data_full$alpha_a, lty = "dashed")
plot(Alter, beta_data, type = "l", ylab = "beta")
lines(Alter, est_data_full$beta_a, lty = "dashed")
plot(min(complex_period_data[,1]):max(complex_period_data[,1]), gamma_data, type = "l",
xlab = "Zeitraum", ylab = "gamma", ylim = c(-34, 34))
lines(min(complex_period_data[,1]):max(complex_period_data[,1]),
gamma_vec_full, lty = "dashed")
lines(min(complex_period_data[,1]):max(complex_period_data[,1]), rep(0,41), lty = "dotdash")
plot(min(complex_period_data[,1]):max(complex_period_data[,1]),
gamma_sprung, type = "l", xlab = "Zeitraum", ylab = "gamma")
erg_gamma_sprung = Lee(Zeitraum_full, Alter, cbind(complex_period_data_sprung[,1:2],
complex_period_data_sprung[,3+1]))
lines(min(complex_period_data[,1]):max(complex_period_data[,1]),
erg_gamma_sprung$gamma, lty = "dashed")
lines(min(complex_period_data[,1]):max(complex_period_data[,1]), rep(0,41), lty = "dotdash")
par(mfrow = c(1,1))
dev.off()
# plot the estimation on the wohle data set versus the prediction
pdf("../../1 Doku/graphics/ErrorLeeSimPredict.pdf", width = 10, height = 8)
par(mfrow = c(2,2))
# start of first half:0
plot(Alter, complex_period_data[(1:96)+96*(1-1),3], type = "l", ylab = "Sterblichkeit",
main = "Jahr 0")
lines(Alter, EstRates[,1], lty = "dashed")
lines(Alter, EstRates_full[,1], lty = "dotted")
# end of first half:20
plot(Alter, complex_period_data[(1:96)+96*(20-1),3], type = "l", ylab = "Sterblichkeit",
main = "Jahr 20")
lines(Alter, EstRates[,21], lty = "dashed")
lines(Alter, EstRates_full[,21], lty = "dotted")
# middle of prediction:30
plot(Alter, complex_period_data[(1:96)+96*(30-1),3], type = "l", ylab = "Sterblichkeit",
main = "Jahr 30")
lines(Alter, EstRates[,31], lty = "dashed")
lines(Alter, EstRates_full[,31], lty = "dotted")
# end of prediction:40
plot(Alter, complex_period_data[(1:96)+96*(40-1),3], type = "l", ylab = "Sterblichkeit",
main = "Jahr 40")
lines(Alter, EstRates[,41], lty = "dashed")
lines(Alter, EstRates_full[,41], lty = "dotted")
par(mfrow = c(1,1))
dev.off()
# plot the estimation on the wohle data set versus the prediction
# log plot
pdf("../../1 Doku/graphics/LogErrorLeeSimPredict.pdf", width = 10, height = 8)
par(mfrow = c(2,2))
# start of first half:0
plot(Alter, log(complex_period_data[(1:96)+96*(1-1),3]), type = "l", ylab = "Sterblichkeit",
main = "Jahr 0")
lines(Alter, log(EstRates[,1]), lty = "dashed")
lines(Alter, log(EstRates_full[,1]), lty = "dotted")
# end of first half:20
plot(Alter, log(complex_period_data[(1:96)+96*(20-1),3]), type = "l", ylab = "Sterblichkeit",
main = "Jahr 20")
lines(Alter, log(EstRates[,21]), lty = "dashed")
lines(Alter, log(EstRates_full[,21]), lty = "dotted")
# middle of prediction:30
plot(Alter, log(complex_period_data[(1:96)+96*(30-1),3]), type = "l", ylab = "Sterblichkeit",
main = "Jahr 30")
lines(Alter, log(EstRates[,31]), lty = "dashed")
lines(Alter, log(EstRates_full[,31]), lty = "dotted")
# end of prediction:40
plot(Alter, log(complex_period_data[(1:96)+96*(40-1),3]), type = "l", ylab = "Sterblichkeit",
main = "Jahr 40")
lines(Alter, log(EstRates[,41]), lty = "dashed")
lines(Alter, log(EstRates_full[,41]), lty = "dotted")
par(mfrow = c(1,1))
dev.off()
return(c(mean(errors), mean(errors_full)))
}
GoodnessLeeSimSprung = function(){
# Estimate the goodness of the sprung data set; is the sprung visible in the estimator?
# estimate parameters and errors
Zeitraum = min(complex_period_data_sprung[,1]):max(complex_period_data_sprung[,1])
Alter = min(complex_period_data_sprung[,2]):max(complex_period_data_sprung[,2])
est_data = rep(NA, length(Alter)*length(Zeitraum))
m = 10000
gammasprung = matrix(rep(NA, m*length(Zeitraum)), nrow = m)
for(i in 1:m){
erg_gamma_sprung = Lee(Zeitraum, Alter, cbind(complex_period_data_sprung[,1:2],
complex_period_data_sprung[,3+i]))
gammasprung[i,] = erg_gamma_sprung$gamma
}
est_gamma_sprung = rep(NA, length(Zeitraum))
for(i in Zeitraum){
est_gamma_sprung[i+1] = mean(gammasprung[,i+1])
}
devtools::use_data(est_gamma_sprung, overwrite = T)
}
KI_nu = function(){
# estimate the goodness of the Lee-Carter Model as predictor; to this end use the first
# half of the simulated data set to predict the second half and compare to the true values
firsthalf = round(min(complex_period_data[,1]) +
length(min(complex_period_data[,1]):max(complex_period_data[,1]))/2)
Zeitraum = min(complex_period_data[,1]):firsthalf
Zeitraum_full = min(complex_period_data[,1]):max(complex_period_data[,1])
Alter = min(complex_period_data[,2]):max(complex_period_data[,2])
# calculate the mean of EstRates
m = 100
EstRatessum_min = matrix(rep(0, length(Alter)*length(Zeitraum_full)), ncol = 41)
EstRatessum_max = matrix(rep(0, length(Alter)*length(Zeitraum_full)), ncol = 41)
EstRatessum = matrix(rep(0, length(Alter)*length(Zeitraum_full)), ncol = 41)
for(i in 1:m){
Deathrates = cbind(subset(complex_period_data[,1:2], complex_period_data[,1] <= firsthalf),
subset(complex_period_data[,3+i], complex_period_data[,1] <= firsthalf))
erg = Lee(Zeitraum, Alter, Deathrates)
ones = matrix(rep(1, 1 + max(complex_period_data[,1]) - min(complex_period_data[,1])),
nrow = 1)
alpha_matrix = matrix(erg$alpha_a, ncol = 1) %*% ones
beta_vec = matrix(erg$beta_a, nrow = 1)
YearsToPredict = max(complex_period_data[,1]) - firsthalf
gamma_t = erg$gamma_t
nu = erg$nu
# Add nu here?
sigma_xi = 2.270456
s = (sigma_xi^2)/(gamma_t[1]-gamma_t[length(gamma_t)])
nu_min = nu + 1.96 * (1/(sqrt(s)*sqrt(41+96)))
nu_max = nu - 1.96 * (1/(sqrt(s)*sqrt(41+96)))
# nu standard
nu_vec = seq(from = nu, to = nu * YearsToPredict, by = nu)
gamma_vec = c(gamma_t, gamma_t[length(gamma_t)] + nu_vec)
beta_gamma = t(beta_vec) %*% t(gamma_vec)
EstRates = exp(alpha_matrix + beta_gamma)
EstRatessum = EstRatessum + EstRates
# min
nu_vec_min = seq(from = nu_min, to = nu_min * YearsToPredict, by = nu_min)
gamma_vec_min = c(gamma_t, gamma_t[length(gamma_t)] + nu_vec_min)
beta_gamma_min = t(beta_vec) %*% t(gamma_vec_min)
EstRates_min = exp(alpha_matrix + beta_gamma_min)
EstRatessum_min = EstRatessum_min + EstRates_min
# max
nu_vec_max = seq(from = nu_max, to = nu_max * YearsToPredict, by = nu_max)
gamma_vec_max = c(gamma_t, gamma_t[length(gamma_t)] + nu_vec_max)
beta_gamma_max = t(beta_vec) %*% t(gamma_vec_max)
EstRates_max = exp(alpha_matrix + beta_gamma_max)
EstRatessum_max = EstRatessum_max + EstRates_max
}
# now take the mean
EstRates_min = EstRatessum_min/m
EstRates_max = EstRatessum_max/m
EstRates = EstRatessum/m
# plot(0:40, gamma_data, type = "l")
# lines(0:40, gamma_vec_max)
# lines(0:40, gamma_vec_min)
# lines(0:40, gamma_vec)
# lines(0:40, rep(0,41))
pdf("../../1 Doku/graphics/KInu.pdf", width = 10, height = 8)
par(mfrow = c(1,2))
plot(Alter, complex_period_data[(1:96)+96*(40-1),3], type = "l", ylab = "Sterblichkeit",
main = "Jahr 40")
lines(Alter, EstRates_min[,41], lty = "dotted")
lines(Alter, EstRates_max[,41], lty = "dotted")
lines(Alter, EstRates[,41], lty = "dashed")
plot(Alter, log(complex_period_data[(1:96)+96*(40-1),3]), type = "l", ylab = "Sterblichkeit",
main = "Jahr 40")
lines(Alter, log(EstRates_min[,41]), lty = "dashed")
lines(Alter, log(EstRates_max[,41]), lty = "dashed")
lines(Alter, log(EstRates[,41]), lty = "dashed")
par(mfrow = c(1,1))
# maximale Breite
max_breite = max(EstRates_min[,41] - EstRates_max[,41])
# returns 0 if the true model lies within the confidence band
return(sum(EstRates_min[,41] < EstRates[,41]) + sum(EstRates_max[,41] > EstRates[,41]))
dev.off()
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.