読書日記

« 2020年2月 | メイン | 2020年4月 »

2020年3月22日 (日)

　Tポイントとかdポイントとか、ああいうロイヤリティ・プログラムの「ポイント」についての資料。仕事の関係で目を通した。論文じゃないけど、読んだものは何でも記録しておこう、ということで...

翁百合(2019) ポイント経済化について －マクロ経済や金融システムへのインプリケーションを探る－. 日本総研リサーチレポート, 2019-010.

　概観：

• NRIによれば2018年度のポイント発行額は9546億円。これは企業の売上とかにポイント適用率と還元率をかけて推計している。いっぽう金融庁は引当金(負債に属する)の合計額の推移を調べている。
• クレカ・家電量販・携帯電話で発行額が多い。伸びているのはプラットフォーマー企業のポイント。
• ポイントの利用率(失効率)はよくわからない。
• ポイント交換のネットワークは複雑化している。必ずしも交換が進んでいるわけではない。交換自体をビジネス化する動きもある(ネットマイルとか)。
• プラットフォーマーはキャッシュレス決済の経済圏をつくりポイント付加で囲い込みを図っている。
• ポイントを投資に転換する動きもある。これについては金融庁が前向きな方針を示している [←へー]

　ミクロにみると、個社発行ポイント(プラットフォーマーでない企業のポイント)とは、企業が消費者の価格弾力性に注目し、将来の値引き分の決済手段をトークンとして渡すもの。いっぽうプラットフォーマーは、還元額の多くを加盟店に負担させる形でトークンを付与する、というちがいがある。

　プラットフォーマーがポイント還元を積極的に活用するとなにがおきるか？　３つにわけてみていく。

　A. 消費の質の変化。先行研究ではふたつの視点がある。
　A1. ポイントがスイッチングコストとして作用し、ロックインが起きる。この観点からいえばポイント交換はスイッチングコストを下げる。
　A2. 行動経済学的視点でいうと、ポイント付与の表現を工夫したり参照点を変えたりしてフレーミング効果を起こす。[このくだりは先行研究列挙に近い。寺地・近(2011行動経済学), 泉谷(2013放送大の紀要), 中川(2015行動経済学), 山本(2011立教の紀要), Liu(2007 J.Mktg)]
　というわけで、ポイントはロックインにより消費者余剰を低下させている可能性がある。また、キャッシュレス化、データ活用による個別消費提案、消費者行動の非合理化(ポイントを貯めるための貯めるための追加消費とか)が起こりうる。

　B. 価格のカスタム化。先行研究ではポイントを非線形プライシングととらえる見方がある。また、企業は景気が悪くなると囲い込んだ顧客に高めの価格をつける傾向があるので、counter cyclicalな効果があるかもといわれている。
　というわけで、物価指数が実態から離れるかも。なお、プラットフォーマー型ポイントは価格決定力が加盟店にないからcounter cyclicalな効果は小さかろう。

　C. ポイントの疑似貨幣化。先行研究では...[中略]。ただし、仮に年5000億のポイントが決済に使われているとしても、電子マネーの決済金額は約5.5兆、それほど大きくない。また、1単位当たりの価値がばらばらだし交換が限定されているし期限があるし還流しない。貨幣の3機能(価値尺度、流通手段、価値貯蔵)をすべて満たさない。今後のプラットフォーマー型ポイントの進展によってはわかんないけど。
　むしろポイントを貯めること自体に楽しみを見出す人も多いわけで、法定通貨にない魅力があるのかもしれない。貨幣の多様性の展望を拓くものともいえる。
　云々。

　... いっちゃなんだけど、もう少し体裁に気を配ればいいのに、と思った。ちゃんと校正するとか、数式を数式フォントで書くとか。こういうのを読む専門家の方は、そのへんあんまし気にしないのかなあ...

　ポイントについての研究を探すとみつかるのはたいてい購買をアウトカムにとった行動経済学的な話なんだけど(ポイント付与と値引きがどう違うかとか)、そうじゃなくて、ポイント獲得行動そのものといえばいいのか、ケインズのいう貨幣愛ならぬポイント愛といえばいいのか、この資料でいうとCの最後の論点に関心があって、理論面での先行研究を探しているんだけど、見当たらない。このレポートでもアンケートを紹介しているだけだった。ううむ。

論文：マーケティング - 読了：翁(2019) マクロ経済に対する「ポイント」の影響

2020年3月19日 (木)

　コロナ騒動で気が休まらない毎日でありますが、感染者数の推移や予測がメディアにさかんに取り上げられていて、ああ、こういう研究ができたら面白かったろうな... などと無責任なことをぼんやり考えていた次第である。
　で、フォルダを整理していたら、「いずれ読む」PDFのフォルダにそっち方面の概説論文がはいっているのをみつけた。全然覚えていないんだけど、仕事の参考になるかと思って深く考えずに保存していたのだと思う。
　せっかくなので、仕事の待ち時間にめくってみた。

西浦博, 稲葉寿 (2006) 感染症流行の予測：感染症数理モデルにおける定量的課題. 統計数理, 54(2), 461-480.

　まずはSIRモデルから。
　感染する可能性のある人を$S$, 感染している人を$I$, 免疫を獲得したか死んだ人を$R$として、
　$\frac{dS(t)}{dt} = -\beta S(t) I(t)$
　$\frac{dI(t)}{dt} = \beta S(t) I(t) - \gamma I(t)$
　$\frac{dR(t)}{dt} = \gamma I(t)$
$\beta$は感染率。$\beta I(t)$を感染力という。$\gamma$は回復率や隔離率で、$\gamma^{-1}$を感染期間が指数分布に従うと仮定した平均値とする。短期的な流行のモデルなので人口動態は無視する。

　総人口を1として、平衡状態$(S(t), I(t), R(t)) = (1,0,0)$において
　$\frac{dI(t)}{dt} = (\beta - \gamma) I(t)$
が成り立つ。感染初期に感染者は$I(t) \approx I(0) \exp{(\beta -\gamma)t}$と指数関数的に増えるわけだ。病気が集団に侵入可能になるのは$\beta - \gamma > 0$のときである。ここから
　$R_0 = \frac{\beta}{\gamma}$
を基本再生産数という。R noughtと発音する。[←完璧な門外漢なので、このくだりを読んだだけでもう元をとったという感じだ。SIRモデルについては読んだことがあったけど、いまメディアでみかける基本再生産数って全然意味がわかってなかった]

　時刻 $t$までに感染を経験した人を $N(t) = S(0) - S(t)$とすると... [残念ながらこの段落は知識不足で理解できなかったので、メモも省略。シグマ代数ってなに？]

　SIRモデルの拡張。
　$I$をさらに、曝露したけどまだ感染性を持たない状態$E$と感染性を持つ状態$I$にわけると
　$\frac{dS(t)}{dt} = -\beta S(t) I(t)$
　$\frac{dE(t)}{dt} = \beta S(t) I(t) - \epsilon E(t)$
　$\frac{dI(t)}{dt} = \epsilon E(t) - \gamma I(t)$
　$\frac{dR(t)}{dt} = \gamma I(t)$
$\epsilon$は「曝露後に感染性を得る率」。これをSEIRモデルという。
　ほかにもモデルが作れる。免疫が獲得できないSIRSモデルとか。

　基本再生産数の意義として、予防接種率の達成目標について考えると...
　効果$\epsilon$のワクチンを接種率$p$で実施したとき、免疫を持たない人の割合は$1-\epsilon p $だから、接触パターンがhomogeneousなら再生産数は$(1-\epsilon p) R_0$。これを1未満にすれば感染症は終息する。$p$について解くと $p > 1 / \epsilon (1-1/R_0)$だ。[←うぉー おもしれええ]
　$R_0$の推定は、さすがに決定論的に考えてると精度とかを考慮できないので、確率論的なモデルを考えて...[めんどくさいからメモは省略するけど、そんなにものすごく難しい話ではなさそう。ふつうはML推定量が閉形式で書き下ろせるらしい。$\beta$の時間変動とか感受性の異質性とかを考慮したモデルをMCMC推定したりするのもあるのだそうだ。ふへえええええ]

　実際には伝播能力は行動の変化とか公衆衛生施策とかで時間変動するので、そういうのを考慮したのを効果的再生産数($R_t$)という。実際の感染ネットワークから求めたり、理論的に求めたりできる。SARSの流行時にWHOが緊急事態のアラートを出した直後に$R_t$が1未満になったという研究があるそうだ[←へえええええ！]

　感染待ち時間・感染性期間の分布を天下りに与えるのではなく、観察に基づき適切な分布を推定するというアプローチもある。
　感染待ち時間は潜伏期間と似ているので(潜伏期間の途中から感染性があることが多いから、ほんとはちょっとちがうんだけど)、調査で調べた潜伏期間の分布で代用したりする(対数正規分布に従うことが多い)。検疫に必要な目安の時間がわかったりする。[途中から難しくなってきたので大幅中略]
　いっぽう感染性期間の分布を知るのは難しい。感染待ち期間と潜伏期間は似ていることが多いけど、感染性期間と症候性期間はちがう。たとえばHIVの場合、感染した直後にウィルスが多く、いったん下がってAIDS発症の前にまた上がる。[...] 実際のウィルス量と疫学的な感染性との評価指標との相関はまだ確実でなく、今後の課題である。[←へえー。個人レベルの機序と集団レベルの現象の関係が意外にはっきりしないってことなのね... ウィルス量の増減も人によっていろいろってことかしらん]

　接触性の異質性を考慮するという方向の拡張もある。[...中略するけど、ひとつのアプローチはそのなかでは接触パターンが均一であるような複数のsubpopulationを考えるという方法で、これをmultitype epidemicモデルというのだそうだ。よくわかんないけど潜在クラスモデルみたいになるのだろうか...]
　最近関心が集まっているのは、接触関数の頻度分布がべき乗分布に従うとき(スケールフリーなとき)で、感染症の制圧は困難を極める。
　接触頻度に関する定量的説明は性感染症から明らかになるだろうと期待されている。接触回数の計数が観察可能だから。[←これ、誰と何回セックスしたかインタビューで正直に答えるだろうから、って意味だろうか...?]

　あれこれ考えているとどんどん複雑になるので、観察データじゃなくて介入研究でパラメータ推定することも多い。[←うぉー面白い！と思ったんだけど、どういう実験になるのか想像がつかなかった...]

　云々、云々。

　...半分くらいはちんぷんかんぷんだったんだけど、仕事と直接関係しないという気楽さもあって、いやー、面白かった。もっと頭がよくて若かったら、こういうのやりたかったなあ... (← 一時の気の迷い)
　ともあれ、研究者の方々のご奮闘をお祈りいたします。

論文：データ解析(2018-) - 読了：西浦・稲葉(2006) 感染症流行の予測

2020年3月11日 (水)

　仕事の都合であれこれ考えていたらわけがわからなくなってしまったので、ちょっとシミュレーションをやってみたら、さらに新たな疑問が生じ... という話を記録しておく。

背景
　ある量的変数の時系列 $Y_t$ があって、それに影響を及ぼしているであろうなんらかの変数の時系列$X_t$があるとする。なんでもいいんだけど、たとえば売上と広告出稿量とか。
　で、$Y_t$に対する$X_t$の効果の大きさをデータから推定したい。こういうこと、よくありますよね。

　話を簡単にするために、分析者は次のことを知っているとしよう。

• $X_t$の$Y_t$に対する効果は即時的、かつ線形である。つまり、
  　$Y_t = \alpha + \beta X_t + V_t$
  とモデル化できる。
• $X_t$は定常である。なんでそんなことを知っているのかわからんが、なにか実質的な知識があるのでしょう。
• $Y_t$の撹乱項$V_t$は一次自己回帰過程 (AR(1)過程)に従う。つまり、
  　$V_t = \rho V_{t-1} + E_t, \ \ E_t \sim N(0, \sigma_E^2)$
  とモデル化できる。これも、まあ、ありそうな仮定ですわね。時系列の自己相関関数を観察した結果、これAR(1)でいいんじゃね？ と思うことはさほど珍しくないと思う。
• 自己回帰係数$\rho$は$0 \leq \rho \leq 1$である。つまり、$V_t$は非負の自己回帰係数を持つ定常AR(1)過程に従う($0 \leq \rho <1$)か、単位根AR(1)過程に従う($\rho = 1$)かのどちらかである。これもそれほど変な仮定ではないと思う。自己回帰はふつう正だし、もし負だったらデータの観察でそれと気がつくだろう。仮に$\rho >1$なら$V_t$はどっかに飛んでいってしまうはずなので、これも除外できる。でも、$\rho <1$か$\rho = 1$かは簡単には判断できない。

問題
　さて、分析者が関心を持っているのは$\beta$である。ここで疑問なのは、$\beta$の推定誤差は$\rho$とどういう風に関連しているのか、という点である。特に次の点について知りたい。

• Q1. $\rho$が大きいとき、つまり結果変数の撹乱項が高い自己相関を持つとき、$\beta$の推定誤差は大きくなるか、それとも小さくなるか。
• Q2. $\rho = 1$であるのにそれと知らず、$\rho < 1$と仮定するモデルを推定したとき、$\beta$の推定誤差は大きくなるか。

　いずれも、私にとってはちょっと切実な疑問だ。
　Q1についていえば... 仕事のなかで時系列データの分析が生じる際、まずは目的変数の時系列を観察して、これからわざわざ説明変数のデータを集めモデリングする労力は報われるかしらん？ と算段することが多いと思う。そういう場面での手掛かりがほしい。見通しが暗いなら早めに白旗を上げたい。
　Q2についていえば... 理屈の上からいえば、$\rho=1$かそうでないかでAR(1)過程の挙動はがらっとかわる。でも、実データの背後にあるデータ生成プロセスが$\rho=1$かどうかなんて、どうせわかりっこない。だから、きっと$\rho<1$だよねと信じてAR(1)誤差を仮定したり、いや$\rho=1$にちがいないと信じて差分時系列を分析したりするわけである(そんなことないですか？)。でも心には常に不安の影が付きまとう。$\rho$についての仮定をしくじることが、$\beta$の推定にとってどのくらい致命的なのかを知りたい。

方法
　気になって夜も眠れないので(大げさな表現)、簡単なシミュレーションをやってみました。

　$t=-100,\ldots,100$について、データを次のように生成した。
　$x'_t \sim unif(0, 1)$
　$e_t \sim N(0, 1)$
　$v_t = \rho v_{t-1} + e_t$
　$y'_t = x't + v_t$
つまり$\beta = 1$である。で、$x'_t, y'_t$から$t=1$以降を切り出し、それぞれを中心化する。これを$x_t, y_t$とする。
　$\rho$を$0, 0.1, 0.2, \ldots, 1.0$の11通りに動かし、それぞれについて1000個のデータセットを生成した。
　それぞれについてデータセットをひとつ選び$x_t, y_t$を描画すると、こんな感じである。

　モデルの推定方法として、次の6つを試す。

• A. 単にOLS推定
• B. Rのforecast::Arima()でML推定
• C. Rのnlme::gls()でFGLS推定
• D. 状態空間モデルとして定式化し、RのKFASパッケージを使ってカルマンフィルタで推定
• E. Mplusでベイズ推定
• F. Stanでベイズ推定

　推定にあたっては、次の5つのアプローチを試してみる。

• アプローチ1. $\rho = 0$と仮定して、
  　$y_t = \beta x_t + e_t$
  　$e_t \sim N(0, \sigma_e^2)$
• アプローチ2. $|\rho| < 1$と仮定して、
  　$y_t = \beta x_t + v_t$
  　$v_t = \rho v_{t-1} + e_t, \ \ |\rho| < 1$
  　$e_t \sim N(0, \sigma_e^2)$
• アプローチ3. $\rho = 1$と仮定し、$\Delta y_t = y_t - y_{t-1}, \Delta x_t = x_t - x_{t-1}$について、
  　$\Delta y_t = \beta \Delta x_t + e_t$
  　$e_t \sim N(0, \sigma_e^2)$
• アプローチ4. $y_t$について単位根検定を行い、単位根を持たないと判断されたらアプローチ2, 持つと判断されたらアプローチ3に進む。
• アプローチ5: アプローチ2と同じなんだけど、$\rho$の範囲について仮定しない。

　私の乏しい理解では、A(OLS), B(ML), C(FGLS)では撹乱項の定常性が仮定されるので、アプローチ5は選べない。いっぽうD(カルマンフィルタ), E,F(ベイズ)では、$\rho$を明示的に制約しないかぎりアプローチ5になる... というように理解しているのだけれど、ここ、全然自信がない。悲しい。誰か私に教えてくださらないでしょうか。

　まあいいや！とにかくシミュレーションしてみたのであります。

単位根検定
　まず、アプローチ4で必要になる単位根検定の結果について紹介しておく。Dickey-Fuller検定を使った。Rコードはこんな感じ。$y_t$がdfIn$y_centeredにはいっている。

library(urca)

oDFTest <- ur.df(dfIn$y_centered, type = "none", lags = 1)

nD <- if_else(oDFTest@teststat[1,1] < oDFTest@cval[2], 0, 1)

oModel <- lm(y_centered ~ 0 + x_centered, data = dfIn)

library(forecast)

oModel <- Arima(

  y            = dfIn$y_centered,

  include.mean = FALSE,

  order        = c(1, 0, 0),

  xreg         = dfIn$x_centered ,

  method       = "ML"

)

library(nlme)

oModel <- gls(

  y_centered ~ 0 + x_centered, 

  corr = corARMA(p=1, q=0), 

  data = dfIn

)

library(KFAS)

mgZ <- array(dim = c(1, 2, nrow(dfIn)))

mgZ[1,1,] <- as.vector(dfIn$x_centered)

mgZ[1,2,] <- 1

oModel <- SSModel(

  dfIn$y_centered ~ -1 + SSMcustom(

    Z     = mgZ,

    T     = matrix(c(1, 0, 0, 1), nrow = 2), # 最後の要素がrho

    R     = matrix(c(1, 0, 0, 1), nrow = 2), 

    Q     = matrix(c(0, 0, 0, 1), nrow = 2), # 最後の要素がsigma_e

    P1    = matrix(c(0,0,0,0), nrow = 2), 

    P1inf = matrix(c(1,0,0,1), nrow = 2)

  ), 

  H = matrix(0)

)

sub_update <- function(par, model){

    # par: (sigma_eの対数, rho)

    # model: 現在のモデル

    model$T[2,2,1] <- par[2]

    model$Q[2,2,1] <- exp(par[1])

    return(model)

}

oFitted <- fitSSM(

    oModel, 

    inits = c(log(var(dfIn$x_centered)/2), 0.5) , 

    updatefn = sub_update, 

    lower = c(-10, -2),

    upper = c(+10, +2),

    method = "L-BFGS-B"

)

oEstimated <- KFS(oFitted$model)

DATA: 

  FILE = est2.dat;

VARIABLE: 

  NAMES = y x; 

  MISSING=.;

ANALYSIS:

  ESTIMATOR = BAYES;

  BITERATIONS = (2000);

MODEL:

  v by (&1);

  v;

  v on v&1; 

  y on v@1 x;

  [y@0];

  y@0;

data {

  int T;

  vector[T] y;

  vector[T] x;

}

parameters {

  real beta;

  real<lower=-1, upper =1> rho; 

  real<lower=0> sigma;

}

model {

  vector[T] v; 

  v = y - beta * x; 

  v[1] ~ normal(0, sigma/sqrt(1-rho^2));

  v[2:T] ~ normal(rho * v[1:(T-1)], sigma);

}