餡子付゛録゛

ソフトウェア開発ツールの便利な使い方を紹介。

行列を使った計算で、SURの感じを掴んでみよう

計量の教科書には載っているものの非線形モデルには使えず、内生性を制御できないので顧みられることの少ないSUR *1 をふと思い出したので、手順を確認してみました。 1. 2つの方程式からなるSURモデル 理屈はGreeneのEconometric Analysisに詳しく載ってい…

相関行列から一定以上の大きさの変数の組み合わせを抽出

100や1000もある説明変数同士の相関係数を見るはだるいと言う話があって、確かに100×100や1000×1000どころか、10×10の行列でも目視で確認すると見落としが出そうです。相関係数0.5以上をリストするコードを書いてみましょう。 # テスト用のデータフレームを…

よくある操作変数法(IV)のための練習データとOLSによる推定バイアス

内生性について上手くツイートでは問題点を説明できないので、念頭に置いているコードを出します。 # # 教科書例の需要と供給のデータを作る # 真のモデル: # S = 2 + 3*p + 4*z + ν # D = 1 - 1*p + μ # S = D # p:価格, S:供給, D:需要, z:気候か何か, νと…

FE-IV練習用データ生成

dとsの式でpが内生変数、zが操作変数ですが、こんなんで。 noi <- 20 # 個体数 t <- 5 # 観測期間 obs <- t*noi i <- c(t(replicate(t, 1:(noi)))) fe <- runif(noi, min=0, max=100) a0 <- c(t(replicate(t, fe))) a1 <- -1 b0 <- 2 b1 <- 3 b2 <- 4 z <- r…

ゲーム理論で考える、じゃんけんの拡張の数値演算

ゲーム理論で考える、じゃんけんの拡張: ニュースの社会科学的な裏側で使ったコードです。通常のゲーム理論の数値演算で扱われる計算手順に習っていないと言うか、フィーリングで描いたので、何か勘違いがあるかも知れないです。 ##########################…

多重共線性を出してみよう

観測数100ぐらいでも、誤差項の分散など次第では多重共線性が推定結果に影響を与える例を作ってみました。 多重共線性が無いケース set.seed(2103) n <- 100 x1 <- 0.3 * runif(n, min=0, max=10) + rnorm(n, sd=1) x2 <- 0.7 * runif(n, min=0, max=10) + r…

拡張Dickey-Fuller検定を横断面方向に拡張する

同一のデータ生成プロセスを踏んでいる複数の系列 y1, y2, y3, … を、それぞれ個別に単位根検定にかけると、定常であるのに標本数が不足して単位根の棄却に失敗しがちになる一方、c(y1, y2, y3, …) と単純に接続して単位根検定した場合、接続面に断続が生じ…

打ち切りリセット付単位根過程のAR(1)での推定

ちょっと特異なデータを、そこそこ正しく推定するための検討用です。結論は、AR(1)にしておけば無問題でした。 ### データ生成 ### genSample <- function(n=250){ ex <- rnorm(n) x <- rep(NA, length = n) x[1] <- ex[1] * 0.2 for (j in 2:n) { # 係数0.9…

オッズを使ってロジスティック回帰してみる

ロジット・モデルと言うと個票に応用する二項選択モデルを最尤法などで推定することを思い浮かべる人が多いと思いますが、都道府県ごとの率など集計データに一般線形回帰をかける事で推定する事もできます。観測数不明の集計データしか取れない場合も多いの…

観察された相対リスクが全て交絡バイアスによるものだとしても、未測定交絡因子の相対リスクはそんなに高くなくても良いかも知れない

「受動喫煙防止法について論点整理①:受動喫煙による健康リスク・死亡者数の推定はどのくらい信用できるか?」と言うブログのエントリーで、受動喫煙の相対リスク1.22倍が未測定交絡因子によって引き起こされたバイアスだとして、未測定交絡因子の相対リスク…

労働市場シグナリング仮説と外部効果の推定値のシミュレーション

労働市場シグナリング仮説が成立していても、計量分析で推定される外部効果は負になるとは限らない事を確認するための、簡単なシミュレーションです。 # # 労働市場シグナリング仮説に基づく賃金計算関数 # edu_h_s: 大卒比率 # return: 大卒賃金, 高卒賃金 …

最適化レベルで結果が変わってしまうC言語のダメなコード

前世紀に書いたgccで-O1と-O3で実行結果が異なるコードです。インライン展開が作用しました。今のバージョンのgccでどうなるかは確認していません。なお、警告は出ると思いますし、処理系によってはコンパイルが通らない事もあると思います。 #include<stdio.h>int *p</stdio.h>…

Rの関数をラッピング

Rは対話型インターフェイスなので、長々と関数に引数を指定するのは手間隙です。グローバル・オプションを設定したくなります。しかし、発展途上のライブラリが多いため、グローバル・オプションをつけると上手く動作しなくなったり、そもそもグローバル・オ…

Rの文字コードのデフォルト値の設定

Windows版のRはテキストファイルの入出力にCP932と言うかShiftJISを使うわけですが、LinuxやMacintosh OS XとWindowsを併用して使っている人は、UTF-8をデフォルトにしたいかも知れません。read.tableのfileEncodingオプションをいちいち指定するのは面倒で…

マクロ経済学の動的計画法説明用コードのMatlabからRへの移植

一橋大学の阿部教授が院生向けのレクチャー・ノートをアップロードしていたのを見つけ、拝読して勉強させて頂いていたのですが、説明用のコード(P.16--20)がMatlabだったのでRに移植してみました。プロプライエタリなMatlab嫌いな人に役立つかも知れないの…

Rで機械学習(SVM)

人間の直観的な認識力を手軽に模倣するフレームワークである機械学習がもてはやされるようになってもう10年以上経つと思いますが、世間では第3次人工知能ブームでまだまだ注目されているようです。扱いやすいパッケージが拡充されているのはもちろん『情報処…

結婚市場からの離脱を考慮した夫婦の社会的階層格差

@kabutoyama_taro氏に、「ネット界隈が想起する"上昇婚"のイメージ」は・・・数学的にありえないと言われた*1のですが、ネット界隈の仮定に基づけばあり得るので、シミュレーションをして例示したいと思います。数式を書いて検証すると煩雑になりますし、「…

Rの並列処理(doParallel/foreach)で繰り返し処理を中断する方法

確認してみたのですが、試行錯誤してみた結果、結論は「ありません」でした*1。仕方が無いので、OpenMP 4.0以前の代替措置*2と同じ事をしてみたいと思います。速度的には改善になるはずです。 library("doParallel")cl <- makeCluster(detectCores()) regist…

プレビュー画面で特定情報を出力しないためのBloggerテンプレートの書き方

Bloggerで外部JavaScriptを読み込ませるとプレビュー画面だけ妙に時間がかかる謎現象が発生するので、その回避のために、プレビュー画面だけそこの部分の出力をしないことにしました。Bloggerの記事のプレビューを爆速にするカスタマイズ | クラウド番外地を…

JavaでOpenMPもどき

気づくとマルチコア化が一般化したため、数値演算をするときにカジュアルにマルチスレッド対応にしたい状況も随分と増えてきました。計算時間が数分の一になる程度ですが、何時間もかかる計算もあります。GoやFortran 2008のように並列化容易なような言語も…

Rでテキストファイルを読み込んで分解

RでCSVファイルを上手く読み込めないときは、read.tableにあわせてテキスト・ファイルを書き換えた方が手軽な気がしますが、それが都合で出来ないときは、テキスト・ファイルとしてCSVファイルを読み込んでから、カンマで分離することもできます。 con <- fi…

お気軽にAR(1)の構造転換点を区間推定してみる

学術的な意義は無いのですが、時系列データの構造転換点を求めたいときがあります。記述的に転換点を仮定して推定することが多いわけですが、主観的になりやすく言い合いになることがあるからです。統計学者に殺されそうな荒業ですが、尤度を使って構造転換…

特定月の分散が大きいときの大数の法則

別のブログで使った大数の法則を示すシミュレーションのソースコードです。 # 乱数から分析データを作る set.seed(20150209) x.sd <- 16.96534 # 想定標準偏差 n <- 24 # 分析年数 x <- as.numeric(n*12) # 分析データ for(i in seq(1, (n-1)*12, 12)){ x[i:…

Rで各集団における大きい方から2番目の値を調べる

ツイッターで見かけた御題なのですが、(学校の課題などで)ありそうなので、置いておきます。 # お試しデータを作る set.seed(20150208) n <- 30 df <- data.frame(type=c("A", "B", "C")[round(runif(n, min=0.5, max=3.5))], value=60 - ((1:n)-10)^2)# …

時系列データの季節調整をしてみよう

時系列データには、季節バイアスが入りがちです。年中行事はもちろんのこと、天候の変化も周期的に発生します。ゴールデンウィークのある5月よりも6月の方が行楽客が少ないとしても、行楽客が減少し出したとは言えないでしょう。そもそも月ごとに日数も異な…

「年/月」の形式を月末日でDate型に変更

月次データで「年/月」と言う形式(e.g. 2013/11、2014/2)は良く見かけると思いますが、Rのread.table関数でデータフレームに読み込むと文字列型*1になってしまいます。日付型でないと、subset関数などで絞り込むときに不便ですね*2。しかし、as.Date関数は…

Excelでもプログラミングができる事を忘れた経済学徒へ

Excel/VBAで大量のワークシートを一括してCSVで保存できます。 とりあえずコードだけ。Excelブック内にマクロを置く/書く方法は「Excel 2010 VBA の基礎知識」(Excel for Mac 2011)などを参照してください。 Sub J_YAMASAKI() For Each mySheet In Workshe…

パネルデータで全員が一斉に説明変数の値を変化させたときのタイムダミーの係数

立命館大学の柴田悠氏への私信なのですが、パネルデータ分析の固定効果モデルに年次などタイムダミーが入っているときに、ある瞬間に個々が同じような行動をとって、説明変数の値を一斉に変化させたときに、それがタイムダミーの係数として観測されない事を…

統計モデルに観測値と観測値の割り算値を入れたときのバイアス

「データ解析のための統計モデリング入門」の著者の久保拓弥氏へのお返事なのですが、別所の「統計モデルに観測値と観測値の割り算値を入れても問題ない」に対して以下のようなコメントをもらったので、具体的な例を作ってみました*1。 …まあ,みなさん割算…

Rの拡張でライフゲームを作ってみる

わざわざC言語で書くほどのものではなく、単に行列の受渡しを確認しただけですが、Rの中で使えるライフゲームを作ってみたので公開します。 以下のソースコードを LifeGame.c と名づけて保存します。 #include <R.h> #include <Rinternals.h>SEXP LifeGame(SEXP m) { SEXP ans, </rinternals.h></r.h>…