Rによる裁量的発生高の計算方法

実証会計学という財務(管理)会計学ファイナンス、応用ミクロ経済学などの様々な研究領域の知見を援用した学際的な研究分野があります。
その中には、発生主義会計によって生じる会計発生高(accruals)に着目して企業の利益マネジメント行動(earnings management activity)を包括的に捉えようとする研究があります。
ここでいう会計発生高とは費用や収益の見越しや繰延べによって計上される経過勘定項目を集約したようなイメージを抱いていただければ大丈夫だと思います。ですので、発生主義会計を採用することで発生する項目を利用した利益調整は粉飾など法令違反による利益操作のことを全く意味しておらず、企業の1つの経営戦略として位置づけられています。本記事では会計発生高を利用した企業の利益マネジメント行動のうち、発生主義会計により非裁量的に生じる部分を除いた裁量的な部分の計算方法をR言語を使用して以下で例示していきたいと思います。

会計発生高とは以下の等式のように利益からキャッシュ・フローを取り除いた部分として定義されています。

会計発生高=利益-キャッシュ・フロー

この会計発生高は企業で調整をした部分(裁量的発生高)と通常のビジネス活動をしているならば自然発生する部分(非裁量的発生高)に分けることができ、以下のようになります。

会計発生高=非裁量的発生高+裁量的発生高

多くの先行研究では裁量的発生高を算出するために非裁量的発生高を産業・年度ごとにOLSで推定し、その目的変数の予測値からの会計発生高を差し引くことで裁量的発生高を計算しています。
この記事ではR言語を利用して非裁量的発生高の推定モデルの1つであるOCF修正ジョーンズ・モデルの推定コードを書いてみたいと思います。

データ

  • 2009年3月から2014年3月までの決算期間における金融業を除く一般事業会社の年次決算データを用意

財務データの読込み

> library(readr)
> fin <- read_csv("C:/Users/admin/data/fq_200903_201403.csv", 
        col_names = TRUE,locale=locale(encoding="SJIS"),
        col_types = cols(YYYYMM = col_date(format = "%Y/%m")))

ここで読み込んだCSVファイルを確認してみます。

> head(fin)

 NCC	         YYYYMM	        MONTH	TA	AR	SALE	NI	OCF	PPE
0000001	         2009-03-01	12	61184	16880	147554	1587	2346	10231
0000002          2009-03-01	12	84026	27504	146273	5092	-3084	16146
0000003	         2009-03-01	12	385462	58515	505250	-16239	-7357	108086
0000004          2009-03-01	12	123511	15214	199239	366	NA	11298
0000006	         2009-03-01	12	11879	2737	35573	272	676	4985
0000015	         2009-03-01	12	52938	1734	54320	1146	2784	38997


各列について説明すると、NCCは企業を識別するコード、YYYYMMは決算年月*1、MONTHは決算月数*2、TAは総資産、ARは売上債権、SALEは売上高、NIは税引後当期純利益、OCFは営業キャッシュフロー、PPEは償却性固定資産です。

年度を表す変数の作成

非裁量的発生高は産業・年度ごとに推定するため、まず年度を表す変数を作成します。そのため、{lubridate}パッケージを使用して以下のようにYYYY変数を作成します。

> library(lubridate)
> fin$YYYY <- year(fin$YYYYMM)
> fin$MM <- month(fin$YYYYMM) #必要はないが月も作成した
> fin <- fin[,c(1:3, 11, 12, 4:10)]
> head(fin,3)

NCC	        YYYYMM	        YYYY	MM	MONTH	TA	AR	SALE	NI	OCF	PPE
0000001         2009-03-01      2009	3	12	61184	16880	147554	1587	2346	10231
0000002         2009-03-01      2009	3	12	84026	27504	146273	5092	-3084	16146
0000003         2009-03-01      2009	3	12	385462	58515	505250	-16239	-7357	108086

産業中分類データの読込み

次に産業を特定するために、ここでは産業中分類コードを使用したいと思います。そのため、企業を識別するコードと産業中分類コードが含まれている新たなデータセットをRに読み込みます。

> nk36 <- read_csv("C:/Users/admin/data/nk36code.csv", 
        col_names = TRUE,locale=locale(encoding="SJIS"))
> head(nk36)

NCC	        NK36
0000001       35
0000002       35
0000003       35
0000004       35
0000005       35
0000006       43

3列目の「NK36」が産業中分類コードです。ここで、データフレームnk36に含まれている企業に重複があるかを確認してみます。具体的にはNCCを主キーとみなして、これが一意なキーとなっているかを確認します。

> head(table(nk36$NCC), 28)

0000001 0000002 0000003 0000004 0000005 0000006 0000011 0000012 0000013 0000014 
      1       1       1       1       1       1       1       1       1       2 
0000015 0000016 0000017 0000018 0000019 0000021 0000022 0000023 0000025 0000026 
      1       1       1       1       2       1       1       1       1       1 
0000027 0000028 0000029 0000037 0000041 0000042 0000043 0000044 
      1       1       1       1       1       1       1       2 

これは、データフレームnk36について、企業を識別するコードNCCごとにデータフレーム内に同一のコードが何件存在するかを表示しています。結果をみると、NCC={0000014, 0000019, 0000044}が2とカウントされており重複していることが分かります。ここでは重複していないものを含めて表示されているので、次に重複しているNCCのみをピックアップしてみます。

> dupli <- nk36$NCC[duplicated(nk36$NCC)]
> head(dupli)
> str(a)

'0000014' '0000019' '0000044' '0000087' '0000128' '0000150'
 chr [1:875] "0000014" "0000019" "0000044" "0000087" "0000128" ...

結果をみると、NCC=={'0000014' '0000019' '0000044'…}とtable関数で2とカウントされているNCCのみピックアップされていることが分かります。これが重複しているNCCの数であり、全部で875個あることが分かります。
さて、この重複しているデータを削除したいと思います。そのため、はじめに現在のnk36の大きさを確認すると、

> dim(nk36)[1]
  6450

となり、ここから重複している875を引いた5575が重複していない企業数となります。なので、以下のようにdistinct関数を使用してnk36を上書きします。

> library(magrittr) #%<>%を使用するため
> nk36%<>% distinct(NCC,.keep_all=TRUE)
> dim(nk36)[1]
  5575

財務データとの結合

次に財務データと産業中分類コードを結合します(NCCを主キーとする)。結合には以下のようにleft_join関数を使用します。

> data <- left_join(fin, nk36 ,by="NCC")
> data <- data[,c(-13)]
> names(data)[1]<-"NAME"
> data <- data[,c(1:2,13,3:12)]
> data %>% select(-1) %>% head()

NCC	NK36	YYYYMM	        YYYY	MM	MONTH	TA	AR	SALE	NI	OCF	PPE
0000001	35	2009-03-01	2009	3	12	61184	16880	147554	1587	2346	10231
0000002	35	2009-03-01	2009	3	12	84026	27504	146273	5092	-3084	16146
0000003	35	2009-03-01	2009	3	12	385462	58515	505250	-16239	-7357	108086
0000004	35	2009-03-01	2009	3	12	123511	15214	199239	366	NA	11298
0000006	43	2009-03-01	2009	3	12	11879	2737	35573	272	676	4985
0000015	71	2009-03-01	2009	3	12	52938	1734	54320	1146	2784	38997

モデルを推定するために必要な変数の作成

非裁量的発生高を推定するために、OCF修正ジョーンズ・モデルを利用します。この推定モデルは1999年のJournal of Accounting Researchに掲載されたKasznikの論文で利用されたものです。詳細についてはググっていただくとして、ここではざっくりと説明を進めます。まず、推定モデルの目的変数と説明変数を作成し、モデルを産業・年度のグループごとに推定するときに利用するグルーピング変数(intra変数)を以下のとおり作成します。

library(dplyr)
samp <- data
# OCF修正ジョーンズ・モデルの説明変数(前期末総資産でデフレート)を作成
> samp %<>% dplyr::group_by(NCC) %>% arrange(YYYY) %>% 
       dplyr::mutate(DSALE = (SALE-lag(SALE)),
             DAR = (AR-lag(AR)),
             d_sar = (DSALE - DAR) / lag(TA),
             ppe = PPE / lag(TA),
             tac = (NI - OCF) / lag(TA),
             d_ocf = (OCF - lag(OCF)) / lag(TA)
  ) 
#産業・年度ごとのグルーピングを表すintra変数を作成
> samp%<>% dplyr::mutate(intra = paste(NK36,YYYY,sep="")) 
> samp$intra <- as.numeric(factor(test$intra)) 

> samp %>% arrange(NCC, YYYY) %>% select(NCC, NK36, YYYY, intra, tac, d_sar, ppe, d_ocf) %>% head()

NCC	NK36	YYYY	intra	tac	        d_sar	        ppe	        d_ocf
0000001	35	2009	103	NA	        NA	        NA	        NA
0000001	35	2010	104	-0.02026674	-0.01211101	0.1890690	-0.0003268828
0000001	35	2011	105	0.07988989	0.21638855	0.1917700	-0.1151615060
0000001	35	2012	106	0.07435814	0.18610335	0.1504582	-0.0028339292
0000001	35	2013	107	-0.07199454	-0.03174117	0.1328161	0.1492988921
0000001	35	2014	108	0.01270947	0.27061085	0.1272989	-0.0657577032

リストワイズ法により推定に用いるサンプルを揃える

推定するにあたり、欠損値について何かしらの処理を行う必要があります。ここでは単純に欠損値のあるサンプルをすべて除去するリストワイズ法を使用します。
(他にもペアワイズ法、平均値代入法、確率的回帰代入法、完全情報最尤推定法、多重代入法など様々な欠損値の補定方法があります)

> wise <- na.omit(samp) #リストワイズ法による処理を施したデータフレームを新たにwiseとする
> dim(wise)[1] # 推定に利用するサンプルサイズ
16813

# 推定に利用する各変数の基本統計量
> wise %>% select(NCC,YYYY, intra, tac, d_sar, ppe, d_ocf) %>% summary()

    NCC                 YYYY          intra            tac           
 Length:16343       Min.   :2010   Min.   :  2.0   Min.   :-6.709265  
 Class :character   1st Qu.:2011   1st Qu.: 63.0   1st Qu.:-0.067786  
 Mode  :character   Median :2012   Median :119.0   Median :-0.037034  
                    Mean   :2012   Mean   :109.7   Mean   :-0.038418  
                    3rd Qu.:2013   3rd Qu.:155.0   3rd Qu.:-0.008089  
                    Max.   :2014   Max.   :198.0   Max.   : 2.182285  
     d_sar               ppe                d_ocf          
 Min.   :-4.59462   Min.   :  0.00005   Min.   :-4.112338  
 1st Qu.:-0.04914   1st Qu.:  0.13875   1st Qu.:-0.032393  
 Median : 0.01390   Median :  0.26878   Median : 0.002756  
 Mean   : 0.02080   Mean   :  0.30172   Mean   : 0.004535  
 3rd Qu.: 0.08315   3rd Qu.:  0.40890   3rd Qu.: 0.039139  
 Max.   : 9.63259   Max.   :157.76038   Max.   :13.722045  

OCF修正ジョーンズ・モデルの推定(非裁量的発生高の推定)

産業・年度ごとに線形モデルを最小二乗法により繰り返し推定するため、以下のようにintraごとにdo関数を利用します。

> ndac_lm <- wise %>% group_by(intra) %>% 
       do(
           fit_ndac = lm(tac ~ d_sar + ppe + d_ocf, data = .)
       ) 

次に推定結果から裁量的発生高を算出するために、まずは{broom}パッケージのaugment関数を利用してサンプルごとの非裁量的発生高(目的変数の予測値)を確認します。

library(broom)
> pred <- augment(ndac_lm, fit_ndac)
> pred %>% head()

intra	tac	    d_sar	    ppe	    d_ocf	    .fitted	    .se.fit	    .resid	    .hat	    .sigma	    .cooksd	    .std.resid
2	-0.08313327	-0.06963683	0.3619978	0.08757376	-0.08112832	0.003334170	-0.002004951	0.01471774	0.02759499	2.017128e-05	-0.07349476
2	-0.08262803	-0.05957199	0.3078039	0.07391090	-0.07287518	0.003434084	-0.009752844	0.01561304	0.02758124	5.072526e-04	-0.35766896
2	-0.09524074	-0.08138149	0.3021316	0.14757325	-0.09902084	0.005361959	0.003780109	0.03806386	0.02759339	1.945516e-04	0.14023749
2	-0.11218558	-0.21276551	0.4864067	0.04740422	-0.07615276	0.004203001	-0.036032823	0.02338756	0.02739744	1.053763e-02	-1.32669184
2	-0.05069609	-0.12227476	0.4567901	-0.03204623	-0.04494688	0.003660201	-0.005749205	0.01773682	0.02759060	2.011135e-04	-0.21107011
2	-0.08836254	-0.51123571	0.1760187	0.08136506	-0.07335015	0.007994533	-0.015012395	0.08461590	0.02755900	7.532706e-03	-0.57092781

サンプルごとの非裁量的発生高は.fitted列に表示されています。したがって、サンプルごとの会計発生高(tac)から非裁量的発生高(.fitted)を引いた値(これが裁量的発生高)を新たにデータフレームpredの列に追加したいと思います。

> pred$dacc <- (pred$tac-pred$.fitted) #predに裁量的発生高を表す変数daccを追加
> summary(pred$dacc)

    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
-1.05300 -0.02171  0.00217  0.00000  0.02440  1.68800 

サンプルごとのdaccも確認してみます。

> pred %>% select(intra, tac, .fitted, dacc) %>% head(5)

intra	tac	        .fitted	        dacc
2	-0.08313327	-0.08112832	-0.002004951
2	-0.08262803	-0.07287518	-0.009752844
2	-0.09524074	-0.09902084	0.003780109
2	-0.11218558	-0.07615276	-0.036032823
2	-0.05069609	-0.04494688	-0.005749205

ここで算出されたdaccや.fittedをもとのデータフレーム(例えばwise)の新たな列に追加します。

> pred %<>% select(intra, tac, d_sar,ppe,d_ocf, .fitted, dacc)
> esti_join <- left_join(wise, pred) # Joining, by = c("intra", "d_sar", "ppe", "tac", "d_ocf")
> esti_join %<>% select(NAME, NCC, YYYY, MM, intra, NI, OCF, tac, .fitted, dacc) %>% rename(ndac = .fitted)
> esti_join %>% select(-NAME) %>% head()

NCC	YYYY	MM	intra	NI	OCF	tac	        ndac	        dacc
0000001	2010	3	104	1086	2326	-0.02026674	-0.03375986	0.013493119
0000002	2010	3	104	2271	12256	-0.11883227	-0.10153250	-0.017299765
0000003	2010	3	104	44	33550	-0.08692426	-0.08024806	-0.006676205
0000006	2010	3	122	227	1370	-0.09622022	-0.06095252	-0.035267696
0000015	2010	3	194	486	1360	-0.01650988	-0.05131384	0.034803958
0000016	2010	3	122	606	2002	-0.04263116	-0.07219979	0.029568625

8列目のndacは非裁量的発生高を表しており、9列目のdaccは裁量的発生高を表しています。
細かな説明はかなり省略して記事を書きましたが、R言語でも非裁量的発生高を推定し裁量的発生高を計算することができました。

*1:もとのcsvでは年月表示だったのでYYYYMMという名前にしたが、読込んでみると年月日となっている。そのためYYYYMMDDとすべきだが今回はYYYYMM表記で記事を書きます。

*2:決算月数が12カ月のデータのみを使用しています。具体的には以下のコードを書いて絞込みました。fin %<>% dplyr::filter(MONTH==12)