R语言：企业风险分析（2）【蒙特卡罗模拟，Monte-Carlo Simulation】

lzc1009840152

已于 2022-04-16 11:21:18 修改

阅读量2.2k

点赞数 1

分类专栏：企业风险分析文章标签： r语言概率论

于 2022-04-11 11:09:16 首次发布

本文链接：https://blog.csdn.net/lzc1009840152/article/details/124074857

版权

企业风险分析专栏收录该内容

5 篇文章 2 订阅

订阅专栏

前言：

在企业风险这一章中，蒙特卡罗方法（蒙特卡罗实验）是一类广泛的计算方法，它依赖于重复随机抽样来获得数值结果。基本概念是使用随机性来解决原则上可能是确定性的问题。Monte Carlo 方法主要用于三个问题类别：优化、数值积分和从概率分布生成绘图。本章主要介绍蒙特卡罗方法的基本操作，计算等。

简介：

一家医药企业正在市场上测试一种新产品：

新产品的需求量为正态分布，平均值为2000000，标准差为250000。预计需求将以每年4%的速度增长。
研发成本估计在5亿至8亿美元之间，最有可能的价值为7亿美元。
临床试验费用估计在1.35亿美元到1.6亿美元之间，最可能的价值为1.5亿美元。
市场上有竞争对手，该企业估计第一年的市场份额将是6%到10%之间的任何数字（可能性相同，该公司估计他们的市场份额将以每年20%的速度增长。
据估计，每月开一次处方可产生240美元的收入。可变成本估计为30美元。

（这里用到了正态分布，三角分布，均匀分布等，有关分布的内容在第一章）

目标：

请开发一个模拟模型，假设贴现率为10%，计算该项目三年内的净现值（NPV）。运行1000次迭代的模拟，并回答以下问题。

a、净现值（平均值和标准差）的分布是什么？
b、净现值的第一个四分位数、中位数和第三个四分位数是多少？
c、三年内净现值不为正的概率是多少？
d、在概率至少为0.9的情况下，我们可能观察到什么NPV？
e、第三年的累积净利润有多大，概率至少为0.9？
f、净现值的95%置信区间是多少？

开始

1. 做好准备工作（创建npv函数），输入已经确定的参数，年份，

library(EnvStats)   #rtri
library(tidyverse)   #%>%

npv <- function(cf,r){
  sum(cf*r^-(1:length(cf)))
}


#10000 trials, 5 years
options(digits=3)
n <- 1000          # 实验次数
y <- 3             # 3年

mkt_size_g <- 1.04   # 市场规模增长
mkt_shr_g <- 1.20    # 市场分额增长

rr <- 1.10             # 利率
u_pft <- (240-40)*12   # 单位利润

2. 创建主模型：

stochastic_model <-
  data.frame(mkt_size <- rnorm(n, mean = 2, sd = .25),
             RD <- rtri(n, 500, 800, 700),
             CT <- rtri(n, 135, 160, 150),
             mkt_share <- runif(n, .06, .10))

效果如图：

3. 创建副表

（因为年份分为1，2，3，对于一次模拟来说，相当于多了一个时间维度，所以创建副表储存和操作）【replicate() 函数常用于创建模拟，它可以重复一个表达式特定次数。】

temp <- replicate(n, simplify = F,
                  expr = data.frame(year = 1:y,
                                    share_growth = c(1,1.20,1.2),
                                    size_growth = c(1,1.04,1.04)))

效果如下：

（从直觉上理解，主表和副表的关系，类似于）

4. 对副表进行操作，计算销售额，每年利润和累计的利润

【cumprod()函数返回累积乘法结果。如cumsum(2:4)，返回结果 ‘ 2 6 24 ’】

【cumsum()函数返回累积加法结果。如cumsum(1:3)，返回结果 ‘ 1 3 6 ’】

for(row in 1:n){
  a <- temp [[row]] %>%
    mutate(sales = stochastic_model$mkt_size[row]*stochastic_model$mkt_share[row]*
             cumprod(size_growth*share_growth),
           profit = sales*u_pft,
           cum_profit = cumsum(profit) - stochastic_model$RD[row] - stochastic_model$CT[row])
  temp[[row]] <- a}

此时副表储存了每次模拟的销售额，每年利润和累计的利润的结果

5. 将副表结果放回主表（注：放回时，因为多一个时间维度，因此处于‘折叠’状态）

stochastic_model$year_df <- temp

6. 开始计算计算NPV

【lapply函数返回与对象的长度相同的列表，每个元素都是应用于FUN的相应元素的结果。】

【sapply函数将列表、向量或数据框作为输入，以向量或矩阵形式给出输出。R中的 Sapply 函数与 lapply函数执行相同的工作，但返回的是向量。】

此处利用lapply打开折叠，得到Profit，并接着sapply得到npv，并减去成本。

stochastic_model$NPV<-
  lapply(X = stochastic_model$year_df,
         FUN = function(df) df$Profit) %>%
  sapply(FUN = npv,r = rr) - stochastic_model$RD - stochastic_model$CT

最后可以用unnest取消year_df的嵌套，把整个表展开

flat_df <- stochastic_model %>% unnest(year_df)

可视化一下：

ggplot()+geom_histogram(aes(x = stochastic_model$NPV, y = ..count.. /sum(..count..)) ,
                        color="dark blue" ,fill="blue" , binwidth=10)

开始答题：

a、净现值（平均值和标准差）的分布是什么？

mean(stochastic_model$NPV)
sd(stochastic_model$NPV)

## [1] -813.4
## [1] 62.64

b、净现值的第一个四分位数、中位数和第三个四分位数是多少？

fivenum(stochastic_model$NPV)

## [1] -935.6 -861.5 -819.5 -767.4 -653.1

c、三年内净现值不为正的概率是多少？

mean((flat_df %>% filter(year==3))$NPV<0)

## [1] 1
# 很不幸，全是负的

d、在概率至少为0.9的情况下，我们可能观察到什么NPV？（左右各5%）

画图示意一下

quantile(stochastic_model$NPV, probs = c(0.05, 0.95))

##  5%  95% 
## -907 -698 

ggplot()+geom_histogram(aes(x = stochastic_model$NPV, y = ..count.. /sum(..count..)) ,
                        color="dark blue" ,fill="blue" , binwidth=10)+
  geom_vline(xintercept = quantile(stochastic_model$NPV, probs = c(0.05, 0.95)),
             linetype='dashed',
             color='red')+
  theme_classic()

e、在概率至少为0.9的情况下，第三年的累积净利润有多大？（与上一题相似）

quantile(filter(flat_df, year==3)$cum_profit, c(.05,.95))

##  5%  95% 
## 182 1150

f、净现值的95%置信区间是多少？

标准方法：

sample.mean <- mean(stochastic_model$NPV)
sample.n <- length(stochastic_model$NPV)
sample.sd <- sd(stochastic_model$NPV)
sample.se <- sample.sd/sqrt(sample.n)
alpha = 0.05       # α是显著性水平
degrees.freedom = sample.n - 1
t.score = qt(p=alpha/2, df=degrees.freedom,lower.tail=F)
margin.error <- t.score * sample.se
lower.bound <- sample.mean - margin.error
upper.bound <- sample.mean + margin.error


lower.bound
## [1] -817
upper.bound
## [1] -810

简便方法：

# Calculate the mean and standard error
model <- lm(NPV ~ 1, stochastic_model)

# Calculate the confidence interval
confint(model, level = 0.95)

##              2.5 % 97.5 %
## (Intercept)  -817   -810

以上为基本的操作和分析，可视化等操作将在下一章（第三章）介绍

内容源于波士顿大学，感谢教授David Ritt

Reference:

apply(), lapply(), sapply(), tapply() Function in R with Examples

Practice 9 Calculating Confidence Intervals in R | R Practices for Learning Statistics