朝が苦手な院生ブログ

基礎こそ物の上手なれ. 人間万事塞翁が馬. を大切にしている経済学徒.

悟りを開く時間が増えてきた気がする(修論提出まであとn週間)

先日、こちらのブログのとある記事がTwitterでプチバズリ(計400-600RTくらい)していたお陰で臨時収入が増えていた。はてなプロを辞めてからアクセス数が激減していたが、全盛期のそれを遥かに超えるアクセス数で、電気代水道代は軽く払える金額なので非常にありがたい。今日は疲れているので、後日そのクソ適当に書いていた記事を書き直そうと思う。
 
PhDの出願は予想を遥かに超える厳しさである。体力的にも、実力的にも。全てを辿れば自分の計画性のなさからくるものである。
 
純粋に修論と並行して行うのは厳しい。M2の冬でPhD出願する人の多くはM1が終わると同時に英語試験から取り掛かるのだと思うが、自分は夏から始めたのでそもそも同じ土俵に立とうとするのが簡単ではない。もちろん競争相手は日本人だけではないが。そもそも、早くから準備していても学部の成績も雑魚だししんどさは変わらないかも。
 
日々ありえないプレッシャーを感じているわけだが、これは多分アプライしてから通知が来るまでの期間の方がしんどいらしい。この一年は、これまでの人生にはなかった挑戦が非常に多く、同じくらい重要なことを学んだ気がする。稚拙な文章だが、色々感じたことを年末に記事にしてみようと思う。私の出願はまだ終わっていないので、もちろんそんなことをしている暇もないのだが、息抜き&一年の振り返りとしての執筆と、高校の部活の友人達との飲み会とスターウォーズ鑑賞だけは心の底から楽しみたい。
 
みなさん、2019年も残り1週間ですが、頑張りましょう。

econgrad.hatenablog.com

「もう疲れた誰か助けてよ!」そんな合図出したって

数週間前にGREのスコアがかえってきた。
 
スコアはVerbalが139点、Quantitativeが162点、AWが4.0点だった。
AW以外カスすぎて笑えない(震え)。
 
テスト1週間前にちょろっと対策したAWは非ネイティブ出願者としてはまあokレベルだと思うのだが、かなり単語勉強したVはゴミクズ、QはEcon出願クラスタの中で言えば下の中という感じだろうか。
クリスマスにもう一度受験するので、最低限Vを140に乗せ、Qを167-168は取りたいところだ。しかし正直なところ、Vの勉強に使える時間など無い。同時に出願に必要なwriting sampleの修論を書きまくっている。僕が在籍する大学が無料で提供しているネイティブの先生によるマンツーマン英語添削サービスを毎週受けにいき、毎回2時間英語の文法ミスやおかしい表現を直して貰っている。しかし、当然ではあるが、「文法的に正しく書けている=いい文章」という訳では無い。指導教授に少しずつライティングの添削も受けている。英語の先生による添削よりも紙が真っ赤になる。
 
他にも某案件(アカハラ)絡みでメンタルが削られて、日々ボロボロになりながらも進んでいる。
 
(この記事の執筆時間は10分です)

修士論文提出まで残り3ヶ月を切った人間

2ヶ月ぶりくらいの投稿です。毎日登校はしているけどね(へっ)。

 
まずこのブログで使っていた独自ドメインの契約をやめました。アドセンスで数千円は稼いだのですが、そもそも頻繁に更新するブログでもなく、良くてプラマイゼロくらいなら必要ないかなと。
 
 

最近は修論の推計結果はほぼまとまり、ライティングを見てもらっているところだ。率直な感想として、楽しいがかなり難しい、そして疲れる。アメリカ某大学の教授が私の修論と同じテーマのWPを2017年から公開していたのですが、それがつい数週間前に某ジャーナルにアクセプトされていた。そもそも経済学は他の学問と違い、二番煎じにとことん厳しい(評価されない)分野だと思うのだが、この論文が公開されてしまったので三番煎じになってしまった点がショックである。開発経済学だと、「この国ではこの研究がされていますが、この国では初めてです。だから重要です。」 みたいな。  
しかし、修士論文はそこまでのliteratureへの貢献は求められないので、まあ残念だがいい経験になったかなと。投稿するジャーナルのランクを下げて、新しい研究に移るのが最適戦略じゃないかなと言われた。修論としてはそこそこユニークなデータセットを構築し(死ぬほど時間がかかった)、robustな推計結果およびfalsification testなどrobustness checksをクリアすることが出来たので、causal inferenceとまでは言えないものの、ミクロ計量の基本的なスキルを身に付けることができたのが大きな成長なのかなと思う。
 
フィールドトップに数本出している先生の指導のもと書いているのだが、ライティングは本当に奥が深いなと感じる。イントロは論文で一番難しく、かつ最初の数行で読者の関心を引かなければいけない。ここで失敗すると、途中で飽きられて、著者がいい研究者ではない可能性があると思われる。先生は昔、イントロを書くのに数ヶ月費やしていたそうだ。何十回もイントロを書き直すうちに、いいものが書けるようになるのだと。そのような努力の積み重ねがトップ誌への掲載に繋がるのかと実感することができた。
 
 
4月から博士課程に在籍する予定なのだが、同時にPhD留学の準備も進めている。進めているつもりだが、やはり準備を始めるのが遅かったこと、修論との同時並行は少し厳しいと感じることがある。あと根本的な問題として、成績(GPA)が微妙であること、推薦状が弱いことが問題だと言われた。これは自分がどのレベルの大学(プログラム)を目指すかによって変わるのだが。研究室の先輩が米トップ30-50くらいのプログラムに進学しているのだが、その先輩と比べても自分の出願書類は弱いが明らかなのが辛い。他人とは比較せず、準備に手を抜かないことが必須だなと改めて感じた。
 
 
疲れたので息抜き程度に20分で日記程度につけてみた。近いうちに、これまでの留学準備を通じて感じたこととかを書こうと思う。

順調な人生と「欲望と必要」という存在

最近こんなことしているよ

気が付けば貴重な貴重な夏休みの半分が終わってしまっているではないか。半年ちょい前からいつ修理に出すか迷っていたMacBookが、先日、問題だった箇所とは違うストレージの問題が起こりQGIS上の計算が回らなくなってしまった。リンゴ社からは既にMacBookの生産はされておらず、MacBook Proへの乗り換えのタイミングを伺っていたので思い切って購入してやった。カスタマイズでメモリを16GBに上げたので値段も+2万ジャポニペリカ払った。あと数日で届く。あとは手元のクソMacBookAppleCareを使って修理に出してから売る。iPad Proも売る。7〜8割はcompensateしたい。クソ金ない。
 
 
最近、人生は驚くほど順調だ。順調に超絶ハードモードで壁が多すぎてヘコたれている。よく聞く話だが、人生の7割8割は辛いらしい。ならツラミが深いことがあると順調なんじゃねえか。と思い始めてきたのである。もしかすると俺は何か悟ったのかもしれない。因みに、ソフトバンクの柳田選手は記者に「悟りましたね。」と言われて「なんスカ悟りましたって?」と返していた。そいうい悟りもあるらしい。
 
 
冒頭で書いたように、MacBookがイかれているためこの1週間研究を放置プレイして楽しんでいる。訳わからない、どうせ点数も上がらんのだろうと早くも捨てモードになりながらも某ぼったくり英語試験に貢いでいる。人生とは惨めだ。惨めすぎてインテリぶりたくなったので、本棚にある本を何冊か読み返した。自分は新しい本を読むのも好きだが、同じくらい本を読み返すのも好きだ。音楽と映画はなかなか新しい作品には手を出しにくい。同じものばかり見ている。コナンとスターウォーズは幼少期から全て見ている。特に『世紀末の魔術師』は神作すぎてエンディングの時点で毎回号泣はしないが全俺が感動の渦に巻き込まれている。以下、半年ちょいぶりに読み返した本について語る。
 
 

毎回のごとく遅めに本題に入るよ

本題に入ろう。今日は「欲望と必要」という存在というものについて語りたい。  
 
「欲望(ウォント)と必要(ニーズ)」というのはフィールズ賞受賞者である数学者・広中平祐著『学問の発見 数学者が語る「考えること・学ぶこと」』から拾ってきた。この本は元々30年以上前に刊行された本で(こちらがオリジナルの本)昨年2018年7月ににブルーバックスから再出版された。
 
広中先生は本書の一部で、人生においてこの「欲望(ウォント)と必要(ニーズ)」の重要性を述べられており、非常に共感したので以下いくつか引用・抜粋しながら感想を述べたい。
 
『「必要」は、英語で主に2通りの表現の仕方がある。ニーズとウォントである。(略)「ニーズ」という言葉は、空間的に言えば、外部の状況を判断して割り出した必要性であり、時間的に見ると、過去から現在にかけて人間が経験したこと、得たものを基準にして割り出した必要性という意味に使われる。これに対して「ウォント」は、自分の内部から出てくる必要性であり、現在と未来に時間軸をとった上での必要性を意味している。(PP161-162)』
 
なるほど。例として、企業の広告などで散見する「お客様のニーズをよく捉えて・・・」などという表現は微妙で、そんなことをしていたらその企業は立ち遅れてしまうという。「お客様のウォントを見抜いて・・・」が正しいのではないか、と言われている。確かに頷いてしまう。
 
自分が強く共感したのは、我々が進路選択など将来の意思決定においてもこの考えを当てはめるべきであるという点だ。
 
『自分の将来を決めていくという時に、いろいろな情報がある。例えば、自分の偏差値がこの程度だからあの大学のこういう学部に行こうとか、こういう職種が有望だからこの企業に就職しようという具合に、いろいろな情報からニーズを割り出して進路を決める人が非常に多い。しかし、そういう決め方をした人は何らかの方法でニーズから割り出したものが、ウォントに切り替わらないかぎり、どこかで挫折するのではないかと思う。「自分はこの学問がしたいんだ」「私はこの仕事につきたいんだ」というウォントを思った意志力がなければならないのである。(PP162-163)』
 
自分も進路にいろいろと迷い、一度留学ではなく就職をするべきだと考えて計2ヶ月ほど就職活動をかじったのだが、それを辞める決め手となったのはこの考えが根底にあった気がする。自分の性格的に好きでないものはとことん出来ないし、この貴重な20代をとりあえず安定しそうだから就職するといった判断で浪費してしまうのが怖くなったのだ。三毒の煩悩の1つである貪欲に駆りたてられているのか、俺はマジョリティとは違う生き方をしたい、することによって人生を楽しめているのである。「成功=高い山に登りきる」といわれるように、人とは違う景色を見ることに満たされる気分になれるのだろう。
 
我々が目標に向けて何か行動をしているとき、頻繁に苦しみを味わい、自己効力感や自己肯定感が保てなくなり、仲間がいないと孤独感に苛まれる。1つの原因は、設定している目標が高すぎて、少ない歩数でそこに到達しようとすることだろう。高い目標設定も最短で辿り着こうとすることも、どちらもよい選択だと思う。しかしそういった目標には一本道で辿り着くことが出来ないだろう。ほぼ100%無理だ。たくさんの紆余曲折を経て、成長だけでなく時には後退しかしないような時期も経験していく必要があると思う。そのような過程の中で、『自分自身のウォントだと思っていたものが、実は、社会の風潮とか、流行とか、あるいはマスコミがもたらす情報とか、そういう形から形成されたウォント』であるとどこかで必ず挫折してしまうのではないか。
 
先に述べたように、一気に高い山に到達するとすぐ息切れしてしまうだろう。ほとんどの場合、自己の努力による成長に不連続な飛躍なんて存在しないと考える。そういったことを踏まえると、このウォントを見つけることは改めて重要なのだと思わされた。もう一度、自分のウォントについても考え直したい。そして遠回りをしながらも反例を積み重ねて、知的に放蕩する人生を歩んでいきたい。と私は悟ったのであった。
 
 
 
そして忘れずにアフィリンクも貼っておく。一人暮らしはウ◯コは大学でするのと同じくらい常識だ。

人生の夏休みの夏休みの過ごし方

大学院にも夏休みというものが存在し、そして私は大学院生であるため、同世代の友人たちは汗水流して働いている中、私は夏休みをしている。
 
 
しかし、そもそも実験等の研究室ワークが無い経済学M2の大学院生には夏休みなどあまり関係ない。自分のペースで自分の修論を進めれば良いのである。学期中もそうだろう。幸い(今のところは)修論は順調に進んでおり、先生からもこれだと「某ジャーナル(フィールドtop5くらい)には載るかもね」と言って貰っているので修論に関してはあまり焦りなどはない。
 
 
それなのに、毎日12-14時間ほど研究室に篭ってヒソヒソと本と向き合っている。
 
 
その本のタイトルにはGREとIELTSと書かれている。そう、お金を払ってテストを受けさせて貰うのだ。
 
 
某「TOEICとか5,000円くらいだよね。そのテスト、聞いたことないけど、有名?いくらするの?」
 
俺「1回の受験で2まん5せんえんくらいするよ。有名だよ、一部界隈では...」

RでStataと同じクラスターロバスト標準誤差を取り扱う

RでStataと同じ様にクラスタロバスト標準誤差を取り扱うには

以前、計量経済学の実証コースで、「実証論文を1本選んで自分でデータを取ってきてReplicateせよ(論文と同じ分析結果を得よ)」という課題があり、この課題で苦労したクラスタロバスト標準誤差の取り扱いについてこの記事を書こうと思う。
 
計量経済学の分析手法を使う研究を行う多くはStataユーザーであると思う(近年、RとPythonユーザーがかなり増えているようだが)。私もB3から計量分析を学び始め、Stata畑で育った人間だ。授業はRで実証分析を学ぶプログラムで、Replicationも当然Rで行い、自分の書いたコードとstargazerなどでまとめた結果を提出せよといったものであった。そこで母語がStataの私が長時間解かす羽目になったの原因が、Rでクラスタロバスト標準誤差を取り扱う方法が分からなかったからである。  

ロバスト標準誤差

Stataのデフォルトでは、誤差項の均一分散を仮定して標準誤差が計算されるので、reg y x1 x2, robustとして標準誤差をheteroskedasticity robustにする。これで不均一分散に対して頑健な標準誤差(heteroskedasticity robust standard error)を考慮した回帰分析を行うことができる。Rで同じ分析を行う方法も少しググれば出てくる。
 

クラスタロバスト標準誤差

Cluster-robust standard error(Cluster-robust ES)をRで推定するにはどうすればいいのだろう?
Cluster-robust standard errorとは、村や地域のクラスター間の誤差項の相関を考慮したSEのことである。例えば、推定された通常の標準誤差が過小推定されていた場合、本来(Cluster-robust ESが正しいとき)は有意でない係数が有意であると判断されることを避けることができる。
 

Stataでは、reg y x1 x2, robustだけではクラスター間の誤差項の相関を無視してしまうため、reg y x1 x2, vce(cluster id)を指定する必要がある。確か、, cluster(cluster id)でも出来た気が...(記憶力ゴミなのでここはご自分でご確認下さい)。
 
これと同じ作業を行うのに少し苦労した。Rのestimatrパッケージを使えばStataで計算されたCluster-robust標準誤差と同じ物を得られる。スペルがestimatorではないことに注意したい。面倒な方法としては、estimatr::coeftest::vcovHCクラスターを指定する方法があるが、出力も中途半端な気がするし、一度推定したものをcoeftest::vcovHCするのが究極的二度手間である(Stataでは一行で書ける)。ここが唯一Rが完全に劣っている点に思えるが、estimatr::lm_robustを使えば一行でクラスタリングするクラスターを指定してCluster-robust ESを計算することが可能になる。
Stataと同じ結果を得る場合には、以下のようにすればいい。

first2_1 <- lm_robust(y ~ x1 + x2, 
                      clusters = cluster id, se_type = "stata", 
                      data = data)

非常に単純なようだが、estimatr::iv_robustで操作変数法を用いる場合には、se_type = "stata"では微妙にStataのそれと異なるSEが計算される(私はそうだった)。
操作変数法の場合には、

first2_1 <- lm_robust(y ~ x1 + x2 | z + x2,
                      clusters = cluster id, se_type = "CR0", 
                      data = data)

と、se_type = "CR0"のように変更すればStataで計算されたCluster-robust ESと同じ結果が得られるはずだ。推定式のzは操作変数、x1が内生変数の場合である。
何故Cluster-robust ESでも通常のlm_robustの場合とiv_robustの場合でse_typeの指定を変更しなければいけないのかは分からない。実際にStataとRの両方を開きながら推定を行ったが、Stataでは論文の分析(Cluster-robust ES)を完全にReplicate出来ているのに、Rではestimatr::iv_robustse_type = "stata"と指定すると、微妙に数値が異なるSEが計算されたのである。
 
おそらく日本語と英語でこの内容(StataとRでの推定結果の違い)が書かれている記事はゼロなので、息の長い記事になることを祈る。次はRmarkdownやLaTexでの資料作成の際に使う参考文献用のBibTexの使い方の記事でも書こうと思う。  
とりあえずReplicateには成功したので安心しているが、何故計算方法(結果)が異なるのか気になる。ご存知の方いらっしゃればコメント欄にて教えて下さい。

社会科学のためのデータ分析入門 章末問題解答(2章-1) Rコード

これまでの章の解答はこちら

章末問題解答(1章-1)と(1章-2)のRコードの記事はこちらこちら(1章-1)(1章-2)から確認できます。
 

はじめに (Textbook Solution: Quantitative Social Science: An Introduction )

Rを使った統計学の日本語のテキストとして非常に定評のある社会科学のためのデータ分析入門の章末問題の解答(Rコード)です。

 
欠点なのかは分かりませんが、こちらのテキストには章末問題の解答がついていません。そして日本語でも英語でもwebで公開されていません(2018年冬ごろの時点では)。2018年冬に私が上巻の章末問題を解いたのですが、一度公開してみようと思ったので複数の記事に分けて投稿していこうと思います。誰かの役に立てればとも思っているのですが、私のコードにミスがあった場合に指摘していただけると嬉しいです。
I would highly appreciate if you could point out mistakes.
 
また同じ変数に関するプロットをする場合でも複数の方法を使ったりもしています。
 

2章-1 (Chapter2 - Section 1)

 
スクリプトをベタ張りしています。

## Chapter 2 Causality 
## Exercise Solution

setwd("~/qss/CAUSALITY") # ご自身のディレクトリを選択

## -----------------------------------------------------
## Taka(the author of this script) uses Japanese-Version QSS.
## -----------------------------------------------------
## Section 1
## Q1

star <- read.csv("STAR.csv")
head(star)

## Create kinder variable.
unique(star$classtype)
star$kinder <- NA 
star$kinder[star$classtype == 1] <- "small" 
star$kinder[star$classtype == 2] <- "middle"
star$kinder[star$classtype == 3] <- "large" 
class(star$kinder)
## Class of tinder has to be factor in the following question.
star$kinder <- as.factor(star$kinder)
unique(star$kinder)

## Re-create race variable.
unique(star$race)
star$race[star$race == 1] <- "white"
star$race[star$race == 2] <- "black"
star$race[star$race == 4] <- "hispanic"
star$race[star$race == 3 | star$race == 5 | star$race == 6] <- "others"


## Q2

small <- subset(star, star$kinder == "small") 
## small <- subset(star, subset = (kinder == "small")) ## same as above
middle <- subset(star, star$kinder == "middle")
large <- subset(star, star$kinder == "large")

## Ignore the missing value by argument na.rm = TRUE
## reading score of 4th grade
sr <- mean(small$g4reading, na.rm = TRUE)
mr <- mean(middle$g4reading, na.rm = TRUE)
lr <- mean(large$g4reading, na.rm = TRUE)

## math score of 4th grade
sm <- mean(small$g4math, na.rm = TRUE)
mm <- mean(middle$g4math, na.rm = TRUE)
lm <- mean(large$g4math, na.rm = TRUE)

## Display & see means of each score.
## Replacing mean by sd, we can get standard deviation of each.
c(sr, mr, lr); c(sm, mm, lm)


## Q3 

srq <- quantile(small$g4reading, 
                probs = seq(0.33, 0.66, 0.33), na.rm = TRUE)
mrq <- quantile(middle$g4reading,
                probs = seq(0.33, 0.66, 0.33), na.rm = TRUE)

smq <- quantile(small$g4math, 
                probs = seq(0.33, 0.66, 0.33), na.rm = TRUE)
mmq <- quantile(middle$g4math,
                probs = seq(0.33, 0.66, 0.33), na.rm = TRUE)

## We can also get them by another way.
## sr33 <- quantile(small$g4reading, probs = 0.33, na.rm = TRUE)
## sr66 <- quantile(small$g4reading, probs = 0.66, na.rm = TRUE)
## 1/3 instead of 0.33 is more acculate.


## Q4

## Make a contingency table.
table(class_size = star$kinder, year = star$yearssmall)

## Function tapply() applies one function repeatedly
## to each level of the factor variable. 
## tapply(x2, x1, mean) means to calculate mean of x2 for x1.

tapply(star$g4reading, star$yearssmall, mean, na.rm = TRUE)
tapply(star$g4reading, star$yearssmall, median, na.rm = TRUE)

tapply(star$g4math, star$yearssmall, mean, na.rm = TRUE)
tapply(star$g4math, star$yearssmall, median, na.rm = TRUE)


## Q5

## White people have higher score in both class size.
tapply(middle$g4reading, middle$race, mean, na.rm = TRUE)
tapply(middle$g4math, middle$race, mean, na.rm = TRUE)

tapply(small$g4reading, small$race, mean, na.rm = TRUE)
tapply(small$g4math, small$race, mean, na.rm = TRUE)


## Q6

tapply(star$hsgrad, star$kinder, mean, na.rm = TRUE)
tapply(star$hsgrad, star$yearssmall, mean, na.rm = TRUE)
tapply(star$hsgrad, star$race, mean, na.rm = TRUE)

 
章末問題解答(1章-1) Rコード
www.econ-stat-grad.com
章末問題解答(1章-2) Rコード www.econ-stat-grad.com

社会科学のためのデータ分析入門 章末問題解答(1章-2) Rコード

章末問題解答(1章-1) Rコードの記事はこちらこちらから確認できます。
 

はじめに (Textbook Solution: Quantitative Social Science: An Introduction )

Rを使った統計学の日本語のテキストとして非常に定評のある社会科学のためのデータ分析入門の章末問題の解答(Rコード)です。

 
欠点なのかは分かりませんが、こちらのテキストには章末問題の解答がついていません。そして日本語でも英語でもwebで公開されていません(2018年冬ごろの時点では)。2018年冬に私が上巻の章末問題を解いたのですが、一度公開してみようと思ったので複数の記事に分けて投稿していこうと思います。誰かの役に立てればとも思っているのですが、私のコードにミスがあった場合に指摘していただけると嬉しいです。
I would highly appreciate if you could point out mistakes.
 
また同じ変数に関するプロットをする場合でも複数の方法を使ったりもしています。
 

1章-2 (Chapter1 - Section 2)

 
スクリプトをベタ張りしています。

## Chapter 1 Introduction
## Exercise Solution

## -----------------------------------------------------
## Taka(the author of this script) uses Japanese-Version QSS.
## -----------------------------------------------------
## Section 2
## Q1

kenya <- read.csv("Kenya.csv")
sweden <- read.csv("Sweden.csv")
world <- read.csv("World.csv")

dim(kenya); dim(sweden); dim(world)

## Create sum of person-year.
kenya$py <- kenya$py.men + kenya$py.women
sweden$py <- sweden$py.men + sweden$py.women
world$py <- world$py.men + world$py.women

## We will create CBRs in 3 ways.
## Create Kenya's CBR. 
kc1 <- sum(kenya[1:15, "births"]) / sum(kenya[1:15, "py"])
kc2 <- sum(kenya[16:30, "births"]) / sum(kenya[16:30, "py"])
kCBR <- c(kc1, kc2)

## Create Sweden's CBR. 
sc1 <- sum(sweden[c(1:15), "births"]) / sum(sweden[c(1:15), "py"])
sc2 <- sum(sweden[c(16:30), "births"]) / sum(sweden[c(16:30), "py"])
sCBR <- c(sc1, sc2)

## Create World's CBR.
## First, we create & add CBR in the data.frame.
world$crb <- world$births / world$py

## Divide data into two periods.  
wp1 <- world[world$period=="1950-1955", ] # first half
wp2 <- world[world$period=="2005-2010", ] # second half

## Calculate CBRs
wc1 <- sum(wp1$births) / sum(wp1$py)
wc2 <- sum(wp2$births) / sum(wp2$py)

## Save CBR as a new vector which has length=2
wCBR <- c(wc1, wc2)

## Print 3 countries' CBRs
kCBR; sCBR; wCBR


## Q2

## We will create ASFR in the different ways again.
## childbearing age: CBA
## Kenya
## 1950-55
Kcba1 <- kenya[4:10, ]
K.ASFR1 <- Kcba1$births / Kcba1$py.women

## 2005-10
Kcba2 <- kenya[19:25, ]
K.ASFR2 <- Kcba2$births / Kcba2$py.women

## total
## K.ASFR <- c(K.ASFR1, K.ASFR2) # easiest way
Kcba3 <- rbind(kenya[4:10, ], kenya[19:25, ])
K.ASFR <- Kcba3$births / Kcba3$py.women
 
## Sweden
## 1950-55
Scbab1 <- sweden[4:10, "births"]
Scbaw1 <- sweden[4:10, "py.women"]
S.ASFR1 <- Scbab1 / Scbaw1

## 2005-10
Scbaw2 <- sweden[19:25, "py.women"]
Scbab2 <- sweden[19:25, "births"]
S.ASFR2 <- Scbab2 / Scbaw2

## total
Scbab <- c(Scbab1, Scbab2)
Scbaw <- c(Scbaw1, Scbaw2)
S.ASFR <- Scbab / Scbaw

## World
## 1950-55
W.ASFR <- world$births / world$py.women
W.ASFR1 <- W.ASFR[4:10]

## 2005-10
W.ASFR2 <- W.ASFR[c(19:25)]

## total
W.ASFR <- W.ASFR[c(4,5,6,7,8,9,10,19,20,21,22,23,24,25)]

## kenya[kenya$age == "15-19", ]

K.ASFR1; S.ASFR1; W.ASFR1
mean(K.ASFR1); mean(S.ASFR1); mean(W.ASFR1)

K.ASFR2; S.ASFR2; W.ASFR2
mean(K.ASFR2); mean(S.ASFR2); mean(W.ASFR2)

K.ASFR; S.ASFR; W.ASFR
mean(K.ASFR); mean(S.ASFR); mean(W.ASFR)


## Q3

## Kenya's TFR
K.TFR <- c(sum(K.ASFR1 * 5), sum(K.ASFR2 * 5))
## the number of women
kw1 <- sum(kenya[kenya$period=="1950-1955", "py.women"]) # first half
kw2 <- sum(kenya[kenya$period=="2005-2010", "py.women"]) # second half

## Sweden's TFR
S.TFR <- c(sum(S.ASFR1 * 5), sum(S.ASFR2 * 5))
## the number of women
sw1 <- sum(sweden[sweden$period=="1950-1955", "py.women"]) # first half
sw2 <- sum(sweden[sweden$period=="2005-2010", "py.women"]) # second half

## World's TFR
W.TFR <- c(sum(W.ASFR1 * 5), sum(W.ASFR2 * 5))
## the number of women
ww1 <- sum(world[world$period=="1950-1955", "py.women"]) # first half
ww2 <- sum(world[world$period=="2005-2010", "py.women"]) # second half
## the number of births
wb1 <- sum(world[world$period=="1950-1955", "births"]) # first half
wb2 <- sum(world[world$period=="2005-2010", "births"]) # second half

## print TFR & #women
K.TFR; S.TFR; W.TFR
kw1; kw2; sw1; sw2; ww1; ww2
wb1; wb2


## Q4

## Kenya's CDR
K.CDR1 <-  sum(kenya[kenya$period=="1950-1955", "deaths"]) / 
  sum(kenya[kenya$period=="1950-1955", "py"]) 

K.CDR2 <-  sum(kenya[kenya$period=="2005-2010", "deaths"]) / 
  sum(kenya[kenya$period=="2005-2010", "py"]) 

K.CDR <- c(K.CDR1, K.CDR2)

## Sweden's CDR
S.CDR1 <- sum(sweden[1:15, "deaths"]) / sum(sweden[1:15, "py"])
S.CDR2 <- sum(sweden[16:30, "deaths"]) / sum(sweden[16:30, "py"])

S.CDR <- c(S.CDR1, S.CDR2)

## World's CDR
W.CDR1 <- sum(world[c(1:15), "deaths"]) / sum(world[c(1:15), "py"])
W.CDR2 <- sum(world[c(16:30), "deaths"]) / sum(world[c(16:30), "py"])

W.CDR <- c(W.CDR1, W.CDR2)

## Print CDRs. 
K.CDR; S.CDR; W.CDR


## Q5

K.ASDR <- kenya[16:30, "deaths"] / kenya[16:30, "py"]
S.ASDR <- sweden[16:30, "deaths"] / sweden[16:30, "py"]

options(scipen = 1) # avoid exponential notation
K.ASDR; S.ASDR
S.ASDR - K.ASDR

## In the every age class Kenya's ASDR is higher. 


## Q5

## rate of Sweden's each age class, P
sp <- sweden$py[16:30] / sum(sweden$py[16:30])

## Kenya's counterfactual CDR
cfK.CDR <- sum(K.ASDR * sp) # result

## Compare the two
K.CDR[2]; cfK.CDR # result: 0.01038914 0.02321646
K.CDR[2] - cfK.CDR

## This means if Kenya has Sweden's population distribution, 
## CDR will be 1 percent point higher than factual CDR. 
## Thus, the comparison of factual CDRs b/w Kenya & Sweden
## does not necessarily give a meaningful result. 

 
章末問題解答(1章-1) Rコード
www.econ-stat-grad.com