diff --git a/Code/Eling_and_Schnell_NAAJ_2020.r b/Code/Eling_and_Schnell_NAAJ_2020.r
new file mode 100644
index 0000000..b942b5d
--- /dev/null
+++ b/Code/Eling_and_Schnell_NAAJ_2020.r
@@ -0,0 +1,74 @@
+# ===================================================
+# 1. Import libraries and SAS OpRisk Global dataset
+# ===================================================
+
+library("rlang")
+library("ggplot2")
+library("MASS") # import plotdist
+library("fitdistrplus") #fitdist
+library("lubridate") # transform to ymd
+library("EnvStats") # import eexp
+library("POT") # import fitgpd
+library("extraDistr")
+library("actuar") # fitdist on log-logistic
+library("sn") # import snormFit
+library("fGarch") # snormFit
+library("evmix") #fgammagpd
+
+sas <- read.csv("./dataset/SAS.csv")
+sas$Month...Year.of.Settlement <- ymd(sas$Month...Year.of.Settlement)
+sas <- subset(sas, Month...Year.of.Settlement < ymd("2015-01-01"))
+sas <- subset(sas, Month...Year.of.Settlement > ymd("1995-01-01"))
+sas <- sas[sas$cyber_risk==1,]
+length(sas$Loss.Amount...M.) # Length is 1593
+colnames(sas)
+
+summary(sas$Shareholder.Equity...M.)
+
+summary(sas$Assets...M.)
+
+# ===================================================
+# Operational Cyber Risk
+# ===================================================
+unique(sas$Basel.Business.Line...Level.1)
+
+insurance_sas <- sas[sas$Basel.Business.Line...Level.1=="Insurance",]
+length(insurance_sas$Shareholder.Equity...M.) #Length is 100
+
+# removing missing values
+insurance_sas <- insurance_sas[!is.na(insurance_sas$Shareholder.Equity...M.),]
+insurance_sas <- insurance_sas[!is.na(insurance_sas$Revenue...M.),]
+insurance_sas <- insurance_sas[!is.na(insurance_sas$Assets...M.),]
+length(insurance_sas$Assets...M.) #Length is 92
+
+# calculate scr_op=min(0.3*bscr, bc)
+insurance_sas$Shareholder.Equity...M.
+insurance_sas$bscr <- 0.3*(insurance_sas$Shareholder.Equity...M./2.11) #0.3bscr
+insurance_sas$bscr
+insurance_sas$bc <- pmax(0.035*insurance_sas$Revenue...M., 0.017*(insurance_sas$Assets...M.- insurance_sas$Shareholder.Equity...M.))
+
+insurance_sas$op <- pmin(insurance_sas$bscr, insurance_sas$bc)
+
+#insurance_sas$op
+#c(min(insurance_sas$op),max(insurance_sas$op), length(insurance_sas$op))
+#c(min(insurance_sas$Loss.Amount...M.),max(insurance_sas$Loss.Amount...M.), length(insurance_sas$Loss.Amount...M.))
+
+plot(insurance_sas$op, insurance_sas$Loss.Amount...M., 
+     xlab = "SCR_op", ylab = "Cyber Losses", 
+     xlim = c(10, 15000),
+     ylim = c(0, 500),
+     main = "Operational Cyber Risk")
+
+# y=x 그래프 추가
+abline(a = 0, b = 1, col = "red")
+
+# 그래프에 각 케이스 점으로 추가
+points(insurance_sas$op, insurance_sas$Loss.Amount...M., col = "black")
+
+cnt <- sum(insurance_sas$op < insurance_sas$Loss.Amount...M.)
+cnt # the number of underestimated cases is 3
+
+# ===================================================
+# 5. Underwriting Cyber Risk
+# ===================================================
+
diff --git a/README.md b/README.md
index af81f9b..e83d049 100644
--- a/README.md
+++ b/README.md
@@ -17,7 +17,7 @@
 ---
 ### 02. Paper Review
 
-#### [Quick-Brief-Summaries](https://github.com/kwoongbae/risk-management-papers/issues)
+### [Quick-Brief-Summaries](https://github.com/kwoongbae/risk-management-papers/issues)
 
 - **규제모형 관련 연구 정리**
 
@@ -25,7 +25,7 @@
     | ------------------------------------------------------------ | ------------------------------------------------------------ |
     | [Eling and Schnell (2020)](https://www.tandfonline.com/doi/abs/10.1080/10920277.2019.1641416) | 데이터 기반으로 scr 산출하고 규제모형 기반으로 scr 산출하여 사이버리스크의 과소평가됨을 보임. |
 
-#### Deep-Dive-Summaries
+### Deep-Dive-Summaries
 
 - **The Drivers of Cyber risk (*JFS 2022*)** [link](https://www.sciencedirect.com/science/article/abs/pii/S1572308922000171) / [review(ko)](./review/Aldasoro_JRS_2022.pdf)
 
diff --git a/code1/Eling_and_Schnell_NAAJ_2020.r b/code1/Eling_and_Schnell_NAAJ_2020.r
new file mode 100644
index 0000000..b942b5d
--- /dev/null
+++ b/code1/Eling_and_Schnell_NAAJ_2020.r
@@ -0,0 +1,74 @@
+# ===================================================
+# 1. Import libraries and SAS OpRisk Global dataset
+# ===================================================
+
+library("rlang")
+library("ggplot2")
+library("MASS") # import plotdist
+library("fitdistrplus") #fitdist
+library("lubridate") # transform to ymd
+library("EnvStats") # import eexp
+library("POT") # import fitgpd
+library("extraDistr")
+library("actuar") # fitdist on log-logistic
+library("sn") # import snormFit
+library("fGarch") # snormFit
+library("evmix") #fgammagpd
+
+sas <- read.csv("./dataset/SAS.csv")
+sas$Month...Year.of.Settlement <- ymd(sas$Month...Year.of.Settlement)
+sas <- subset(sas, Month...Year.of.Settlement < ymd("2015-01-01"))
+sas <- subset(sas, Month...Year.of.Settlement > ymd("1995-01-01"))
+sas <- sas[sas$cyber_risk==1,]
+length(sas$Loss.Amount...M.) # Length is 1593
+colnames(sas)
+
+summary(sas$Shareholder.Equity...M.)
+
+summary(sas$Assets...M.)
+
+# ===================================================
+# Operational Cyber Risk
+# ===================================================
+unique(sas$Basel.Business.Line...Level.1)
+
+insurance_sas <- sas[sas$Basel.Business.Line...Level.1=="Insurance",]
+length(insurance_sas$Shareholder.Equity...M.) #Length is 100
+
+# removing missing values
+insurance_sas <- insurance_sas[!is.na(insurance_sas$Shareholder.Equity...M.),]
+insurance_sas <- insurance_sas[!is.na(insurance_sas$Revenue...M.),]
+insurance_sas <- insurance_sas[!is.na(insurance_sas$Assets...M.),]
+length(insurance_sas$Assets...M.) #Length is 92
+
+# calculate scr_op=min(0.3*bscr, bc)
+insurance_sas$Shareholder.Equity...M.
+insurance_sas$bscr <- 0.3*(insurance_sas$Shareholder.Equity...M./2.11) #0.3bscr
+insurance_sas$bscr
+insurance_sas$bc <- pmax(0.035*insurance_sas$Revenue...M., 0.017*(insurance_sas$Assets...M.- insurance_sas$Shareholder.Equity...M.))
+
+insurance_sas$op <- pmin(insurance_sas$bscr, insurance_sas$bc)
+
+#insurance_sas$op
+#c(min(insurance_sas$op),max(insurance_sas$op), length(insurance_sas$op))
+#c(min(insurance_sas$Loss.Amount...M.),max(insurance_sas$Loss.Amount...M.), length(insurance_sas$Loss.Amount...M.))
+
+plot(insurance_sas$op, insurance_sas$Loss.Amount...M., 
+     xlab = "SCR_op", ylab = "Cyber Losses", 
+     xlim = c(10, 15000),
+     ylim = c(0, 500),
+     main = "Operational Cyber Risk")
+
+# y=x 그래프 추가
+abline(a = 0, b = 1, col = "red")
+
+# 그래프에 각 케이스 점으로 추가
+points(insurance_sas$op, insurance_sas$Loss.Amount...M., col = "black")
+
+cnt <- sum(insurance_sas$op < insurance_sas$Loss.Amount...M.)
+cnt # the number of underestimated cases is 3
+
+# ===================================================
+# 5. Underwriting Cyber Risk
+# ===================================================
+
diff --git a/code1/Eling_et_al_EJOR_2019.r b/code1/Eling_et_al_EJOR_2019.r
new file mode 100644
index 0000000..3522bac
--- /dev/null
+++ b/code1/Eling_et_al_EJOR_2019.r
@@ -0,0 +1,86 @@
+# ===================================================
+# 1. Import libraries and SAS OpRisk Global dataset
+# ===================================================
+
+library("rlang")
+library("ggplot2")
+library("MASS") # import plotdist
+library("fitdistrplus") #fitdist
+library("lubridate") # transform to ymd
+library("EnvStats") # import eexp
+library("POT") # import fitgpd
+library("extraDistr")
+library("actuar") # fitdist on log-logistic
+library("sn") # import snormFit
+library("fGarch") # snormFit
+
+sas <- read.csv("./dataset/SAS.csv")
+sas$Month...Year.of.Settlement <- ymd(sas$Month...Year.of.Settlement)
+sas <- subset(sas, Month...Year.of.Settlement < ymd("2014-03-01"))
+sas <- subset(sas, Month...Year.of.Settlement > ymd("1995-01-01"))
+sas <- sas[sas$cyber_risk==1,]
+length(sas$Loss.Amount...M.) # Length is 1593 (In Eling 2019, length is 1579)
+
+# ===================================================
+# 2. Goodness-of-fit analysis : severity
+# ===================================================
+
+loss <- sas$Loss.Amount...M.
+plotdist(loss)
+
+# Fit on Exponential distribution
+exp.model <- eexp(loss, method = "mle") # rate = 0.02770837
+exp.loglik <- sum(dexp(x=loss, rate = exp.model$param["rate"], log=T))
+exp.loglik # -7305.53 (In paper, -7535.78)
+
+# Fit on Gamma distribution
+gam.model <- egamma(loss, method = "mle") # shape=0.2416648, scale=149.3397808
+gam.loglik <- sum(dgamma(x=loss, shape = gam.model$param["shape"], scale = gam.model$param["scale"], log=T))
+gam.loglik # -5183.15 (In paper, -5368.23)
+
+# Fit on GPD distribution
+gpd.u <- quantile(loss, probs = 0.9, type = 1)
+gpd.model <- fitgpd(loss, gpd.u, "mle") #scale=60.9383, shape=0.9742
+gpd.loglik <- sum(
+    dgpd(x=loss, gpd.u, gpd.model$param["scale"], gpd.model$param["shape"], log = T)
+)
+gpd.loglik # -Inf (In paper, -4553.42)
+
+# Fit on Log-logistic distribution
+llogis_model <- fitdist(loss, "llogis", method = "mle") #shape=0.868823, scale=1.591586
+llogis.loglik <- llogis_model$loglik
+llogis.loglik # -4399.835 (In paper, -4588.09)
+
+# Fit on Weibull distribution
+wei_model <- fitdist(loss, "weibull", method = "mle")
+weibull.loglik <- wei_model$loglik
+weibull.loglik # -4700.99 (In paper, -4886.78)
+
+# Fit on Skew-normal distribution
+snorm_model <- snormFit(loss)
+snorm.loglik <- sum(
+    dsnorm(x=loss, 
+           mean=snorm_model$par["mean"], 
+           sd=snorm_model$par["sd"],
+           xi=snorm_model$par["xi"], log=T)
+)
+snorm.loglik # -10503.05
+
+# ===================================================
+# 2.1 Goodness-of-fit analysis : severity with POT 56%
+# ===================================================
+
+loss.56 <- sort(loss)[1:round(length(loss)*0.56)]
+loss.44 <- sort(loss)[round(length(loss)*0.56+1):round(length(loss))]
+plotdist(loss.56)
+plotdist(loss.44)
+
+# Fit on GPD distribution (Tail 44% distribution)
+gpd.u <- quantile(loss, probs = 0.59, type = 1)
+gpd.model <- fitgpd(loss.44, gpd.u, "mle") #scale=60.9383, shape=0.9742
+gpd.loglik <- sum(
+    dgpd(x=loss, gpd.u, gpd.model$param["scale"], gpd.model$param["shape"], log = T)
+)
+gpd.loglik # -Inf (In paper, -4553.42)
+
+