R语言：企业风险分析（3）【蒙特卡罗模拟的可视化，敏感性分析】

简介：

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

本章主要介绍蒙特卡罗方法的可视化（蒙特卡罗方法的基本操作，计算在上一章）

开始前基本处理【如若不理解某项操作，请看上一章】

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   # 单位利润

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))

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)))

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}



stochastic_model$year_df <- temp

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

flat_df <- stochastic_model %>% unnest(year_df)

开始相关可视化

1. 叠加图（Overlay Charts）

在一张图表中比较不同年份的累计净利润分布（第1年和第3年），相应地过滤数据，指定分组方式，并填充参数以区分我们正在比较的内容。第3年累计净利润的方差大于第1年累计净利润的方差，这是有道理的，因为进行长期预测时存在更多的不确定性（直观来说，第三年更扁平）

ggplot(flat_df %>% filter(year == 1 | year == 3)) +
  geom_density(aes(x = cum_profit, group = year, fill = as.factor(year)), alpha=0.5)

2. 趋势图（Trend Charts）

创建一个趋势图，比较数据每年的变化情况。首先总结模型，在新的数据框架（summodel）中确定每年75%和90%的平均置信区间。我们把它们命名为l75，u75，l90和u90。每年，我们的累计净利润有75%的几率在l75和u75之间下降，90%的几率在l90和u90之间下降。

summodel <- flat_df %>% group_by(year)%>%   #create a new dataframe 
  summarise(mean = mean(cum_profit),
            l90=quantile(cum_profit, .05),
            l75=quantile(cum_profit, .125),
            u75=quantile(cum_profit,.875),
            u90=quantile(cum_profit,.95))

然后我们在ggplot中绘制它。使用geom_line线和geom_point来表示平均值，geom_ribbon表示两个置信区间。geom_ribbon的alpha参数控制色带的透明度。平均净累积利润随着时间的推移而增加。当在水平轴上向右移动时，75%的频带和90%的频带变得更宽。这表明净累计利润的变化随着时间的推移而增加。

ggplot(summodel)+geom_line(aes(x=year , y=mean)) +
  geom_point(aes(x = year , y = mean)) +
  geom_ribbon(aes(x = year , ymin = l90 , ymax = u90) , 
              alpha=0.3,fill = "blue ")+
  geom_ribbon(aes(x = year , ymin = l75 , ymax = u75),
              alpha=0.3,fill = "blue")+
  theme_bw() +
  ylab("Profit") +
  ggtitle("Trend chart of cumulative Profit") +
  theme(plot.title = element_text(hjust = 0.5))

3. 箱线图(Boxplot charts)

箱线图有助于可视化输出变量的统计特性。在R中，geom_boxplot 将显示第1和第3个四分位数（由方框表示）、中间值（由穿过方框的水平线表示）和垂直线（最多为四分位数范围长度的1.5倍），这些垂直线延伸到方框的上方和下方。任何超过垂直线的数据都被视为异常值。

同样，我们按年份对数据进行分组。

ggplot(flat_df)+geom_boxplot(aes(y=cum_profit, group=year, x=year), outlier.size=0) +
  theme_bw()+
  ggtitle("Boxplot of cumulative Profit") +
  ylab("cumulative profit") +
  theme(plot.title = element_text(hjust = 0.5))

4. 敏感性分析(Sensitivity Analysis)

灵敏度分析的一种形式是测量模型输入和输出之间的相关性。我们将制作一张龙卷风图（tornado chart），用于测量模型中每个变量对净现值的影响。一旦确定了哪些要素具有最高的相关性（正相关性或负相关性），一个明智的组织可能会优先考虑那些对其成功影响最大的投入。

先取出要分析的变量，重命名，得到和NPV的相关性：

c_model <- stochastic_model[-5]
names(c_model) <- c('需求量','研发成本','临床试验费','市场份额','净现值')
cmodel <- cor(c_model)

NPVCorrs <- cmodel[,5][-5]

画图

ggplot() +
  geom_bar(aes(x=reorder(names(NPVCorrs), abs(NPVCorrs)),
               y=NPVCorrs,fill=(NPVCorrs>0)),
           stat= 'identity', color= 'black' )+ 
  coord_flip()+xlab('Random variable')+ 
  ylab('correlation')+
  scale_fill_manual(values=c('red', 'blue'), guide='none')+ 
  theme_classic()+
  ggtitle("correlations between NPv and stochastic elements")+
  theme(plot.title = element_text(hjust = 0.5))
      
# coord_flip :翻转笛卡尔坐标

 显然数据取的不太好，'需求量','研发成本','临床试验费','市场份额'对'净现值'的影响都是负的，因此条带都在左侧。

若一般情况下，图的效果如下：

 

大功告成，快尝试一下吧！

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