R 语言：计算每个公司本年及之前年份的标准差_计算累计标准差

有个小伙伴问了这样一个问题：

如果是求 N 年的投资支出的标准差呢？比如一个公司有 2011～2015 的投资支出数据，2012 年的标准差是用 2011 年和 2012 年的数据，2013 年的标准差是用 2011 2012 2013 的数据，以此类推。

实际上也就是如何计算累计标准差。

我准备了一个示例数据：

library(tidyverse)
haven::read_dta("data1.dta") -> df

df

#> # A tibble: 119 × 3
#>     年份 企业名称                资产总计千元
#>    <dbl> <chr>                          <dbl>
#>  1  2013 ABB新会低压开关有限公司      1185125
#>  2  1998 ABB新会低压开关有限公司        58834
#>  3  2007 ABB新会低压开关有限公司       551123
#>  4  2001 ABB新会低压开关有限公司       121252
#>  5  1999 ABB新会低压开关有限公司        48740
#>  6  2009 ABB新会低压开关有限公司       716935
#>  7  2011 ABB新会低压开关有限公司      1050251
#>  8  2003 ABB新会低压开关有限公司       308216
#>  9  2010 ABB新会低压开关有限公司       917815
#> 10  2000 ABB新会低压开关有限公司        84644
#> # ℹ 109 more rows

按照企业名称和年份排序：

df %>%
  arrange(企业名称, 年份)
#> # A tibble: 119 × 3
#>     年份 企业名称                资产总计千元
#>    <dbl> <chr>                          <dbl>
#>  1  1998 ABB新会低压开关有限公司        58834
#>  2  1999 ABB新会低压开关有限公司        48740
#>  3  2000 ABB新会低压开关有限公司        84644
#>  4  2001 ABB新会低压开关有限公司       121252
#>  5  2002 ABB新会低压开关有限公司       199332
#>  6  2003 ABB新会低压开关有限公司       308216
#>  7  2004 ABB新会低压开关有限公司       342714
#>  8  2005 ABB新会低压开关有限公司       393714
#>  9  2006 ABB新会低压开关有限公司       466189
#> 10  2007 ABB新会低压开关有限公司       551123
#> # ℹ 109 more rows

根据总体方差公式：D(x) = E(x^2) - [E(x)]^2

首先计算 [E(x)]^2：

df %>%
  arrange(企业名称, 年份) %>%
  group_by(企业名称) %>%
  mutate(cumsum = cumsum(资产总计千元),
         n = 1:n(),
         mean2 = (cumsum / n) ^ 2) -> df1

df1
#> # A tibble: 119 × 6
#> # Groups:   企业名称 [7]
#>     年份 企业名称                资产总计千元  cumsum     n        mean2
#>    <dbl> <chr>                          <dbl>   <dbl> <int>        <dbl>
#>  1  1998 ABB新会低压开关有限公司        58834   58834     1  3461439556
#>  2  1999 ABB新会低压开关有限公司        48740  107574     2  2893041369
#>  3  2000 ABB新会低压开关有限公司        84644  192218     3  4105306614.
#>  4  2001 ABB新会低压开关有限公司       121252  313470     4  6141465056.
#>  5  2002 ABB新会低压开关有限公司       199332  512802     5 10518635648.
#>  6  2003 ABB新会低压开关有限公司       308216  821018     6 18724182120.
#>  7  2004 ABB新会低压开关有限公司       342714 1163732     7 27638207507.
#>  8  2005 ABB新会低压开关有限公司       393714 1557446     8 37900594421.
#>  9  2006 ABB新会低压开关有限公司       466189 2023635     9 50556773003.
#> 10  2007 ABB新会低压开关有限公司       551123 2574758    10 66293787586.
#> # ℹ 109 more rows

然后计算 E(x^2)：

df1 %>%
  mutate(x2cumsum = cumsum(资产总计千元^2) / n) -> df2

df2

#> # A tibble: 119 × 7
#> # Groups:   企业名称 [7]
#>     年份 企业名称                资产总计千元  cumsum     n       mean2 x2cumsum
#>    <dbl> <chr>                          <dbl>   <dbl> <int>       <dbl>    <dbl>
#>  1  1998 ABB新会低压开关有限公司        58834   58834     1     3.46e 9  3.46e 9
#>  2  1999 ABB新会低压开关有限公司        48740  107574     2     2.89e 9  2.92e 9
#>  3  2000 ABB新会低压开关有限公司        84644  192218     3     4.11e 9  4.33e 9
#>  4  2001 ABB新会低压开关有限公司       121252  313470     4     6.14e 9  6.93e 9
#>  5  2002 ABB新会低压开关有限公司       199332  512802     5     1.05e10  1.35e10
#>  6  2003 ABB新会低压开关有限公司       308216  821018     6     1.87e10  2.71e10
#>  7  2004 ABB新会低压开关有限公司       342714 1163732     7     2.76e10  4.00e10
#>  8  2005 ABB新会低压开关有限公司       393714 1557446     8     3.79e10  5.44e10
#>  9  2006 ABB新会低压开关有限公司       466189 2023635     9     5.06e10  7.25e10
#> 10  2007 ABB新会低压开关有限公司       551123 2574758    10     6.63e10  9.56e10
#> # ℹ 109 more rows

方差等于两者的差值：

# 方差
df2 %>%
  mutate(var = x2cumsum - mean2) %>%
  ungroup() -> df3

df3

#> # A tibble: 119 × 8
#>     年份 企业名称             资产总计千元 cumsum     n   mean2 x2cumsum     var
#>    <dbl> <chr>                       <dbl>  <dbl> <int>   <dbl>    <dbl>   <dbl>
#>  1  1998 ABB新会低压开关有限…        58834 5.88e4     1 3.46e 9  3.46e 9 0
#>  2  1999 ABB新会低压开关有限…        48740 1.08e5     2 2.89e 9  2.92e 9 2.55e 7
#>  3  2000 ABB新会低压开关有限…        84644 1.92e5     3 4.11e 9  4.33e 9 2.29e 8
#>  4  2001 ABB新会低压开关有限…       121252 3.13e5     4 6.14e 9  6.93e 9 7.84e 8
#>  5  2002 ABB新会低压开关有限…       199332 5.13e5     5 1.05e10  1.35e10 2.97e 9
#>  6  2003 ABB新会低压开关有限…       308216 8.21e5     6 1.87e10  2.71e10 8.35e 9
#>  7  2004 ABB新会低压开关有限…       342714 1.16e6     7 2.76e10  4.00e10 1.23e10
#>  8  2005 ABB新会低压开关有限…       393714 1.56e6     8 3.79e10  5.44e10 1.65e10
#>  9  2006 ABB新会低压开关有限…       466189 2.02e6     9 5.06e10  7.25e10 2.19e10
#> 10  2007 ABB新会低压开关有限…       551123 2.57e6    10 6.63e10  9.56e10 2.93e10
#> # ℹ 109 more rows

不过我们这里实际上是要计算样本标准差，所以需要进行无偏性调整：

回忆一下，总体方差公式里面的分母是 n，样本方差公式里面的分母是 (n - 1)，所以这里需要乘以 n/(n-1) 进行转换：

# 无偏性调整
df3 %>%
  mutate(var = var * (n / (n - 1))) %>%
  select(-cumsum, -n, -mean2, -x2cumsum)
#> # A tibble: 119 × 4
#>     年份 企业名称                资产总计千元          var
#>    <dbl> <chr>                          <dbl>        <dbl>
#>  1  1998 ABB新会低压开关有限公司        58834         NaN
#>  2  1999 ABB新会低压开关有限公司        48740    50944418
#>  3  2000 ABB新会低压开关有限公司        84644   342857025.
#>  4  2001 ABB新会低压开关有限公司       121252  1045940390.
#>  5  2002 ABB新会低压开关有限公司       199332  3710937345.
#>  6  2003 ABB新会低压开关有限公司       308216 10017787511.
#>  7  2004 ABB新会低压开关有限公司       342714 14403243921.
#>  8  2005 ABB新会低压开关有限公司       393714 18813267786.
#>  9  2006 ABB新会低压开关有限公司       466189 24652357071.
#> 10  2007 ABB新会低压开关有限公司       551123 32558722096.
#> # ℹ 109 more rows

点击这里跳转到 RStata 短书平台获取附件：R 语言：计算每个公司本年及之前年份的标准差/计算累计标准差