このページでは、株式会社インサイト・ファクトリー 2020年度社内研修「マーケティング・ミックス・モデリング」VIII章「動学的市場反応モデル(2)」のRコードを公開しています。

準備

library(tidyverse)
## -- Attaching packages --------------------------------------------------------- tidyverse 1.3.0 --
## √ ggplot2 3.3.0     √ purrr   0.3.3
## √ tibble  3.0.0     √ dplyr   0.8.5
## √ tidyr   1.0.2     √ stringr 1.4.0
## √ readr   1.3.1     √ forcats 0.5.0
## -- Conflicts ------------------------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
library(lubridate)
## 
## Attaching package: 'lubridate'
## The following object is masked from 'package:base':
## 
##     date
library(assertr)
library(patchwork)
library(KFAS)
library(dynlm)
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
WORKDIR <- "."

# 「売上と広告」データの作成(Chep6と同じ)
set.seed(123)
dfSalesAds <- tibble(
  # 年月
  dMonth = seq(as.Date("2013-09-01"), as.Date("2018-12-01"), by = "month"),
  # TV広告支出
  gTV    = as.vector(
    stats::filter(
      rnorm(length(dMonth), mean = 5, sd = 2),
      filter = 0.4,
      method = "recursive" # AR
    )
  ),
  # Web広告支出
  gWeb = as.vector(
    stats::filter(
      rnorm(length(dMonth), mean = 2, sd = 1.5),
      filter = 0.3,
      method = "recursive" # AR
    )
  ),
  # 売上の撹乱項
  gYErr = as.vector(
    stats::filter(
      rnorm(length(dMonth), mean = 0, sd = 1),
      filter = 0.45,
      method = "recursive"
    )
  )
) %>%
  mutate(
    # 広告支出の負値を0に置き換え、丸める
    gTV  = round(if_else(gTV < 0, 0, gTV), digits = 1),
    gWeb = round(if_else(gWeb < 0, 0, gWeb), digits = 1),
    # 売上の生成
    gSales = 6 + 0.05 * seq_along(dMonth),
    gSales = gSales + 0.10 * gTV ,
    gSales = gSales + 0.30 * lag(gTV, n = 1) ,
    gSales = gSales + 0.20 * lag(gTV, n = 2),
    gSales = gSales + 0.10 * lag(gTV, n = 3),
    gSales = gSales + 0.05 * lag(gTV, n = 4),
    gSales = gSales + 0.50 * gWeb,
    gSales = gSales + 0.30 * lag(gWeb, n = 1),
    gSales = gSales + gYErr
  ) %>%
  # 変数選択
  dplyr::select(dMonth, gSales, gTV, gWeb) %>%
  # ここで欠損行を削除
  drop_na() 
## print(summary(dfSalesAds))

Code 1. 時間変動つき部分調整モデル

  # データ作成
  tsSales <- ts(dfSalesAds$gSales)
  tsTV    <- ts(dfSalesAds$gTV)
  tsWeb   <- ts(dfSalesAds$gWeb)
  tsDelay <- ts(c(9:11, 0:8)[month(dfSalesAds$dMonth)])

  # 推定
  oModel <- dynlm(tsSales ~ tsTV + tsWeb + tsDelay*tsTV + tsDelay*tsWeb + L(tsSales))
  print(summary(oModel))
## 
## Time series regression with "ts" data:
## Start = 2, End = 60
## 
## Call:
## dynlm(formula = tsSales ~ tsTV + tsWeb + tsDelay * tsTV + tsDelay * 
##     tsWeb + L(tsSales))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.7282 -0.7960  0.0981  0.6144  3.8355 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    4.429928   2.114372   2.095   0.0410 *  
## tsTV           0.167616   0.161391   1.039   0.3038    
## tsWeb          0.408155   0.202337   2.017   0.0489 *  
## tsDelay       -0.019055   0.228771  -0.083   0.9339    
## L(tsSales)     0.594718   0.101640   5.851  3.3e-07 ***
## tsTV:tsDelay  -0.003638   0.024775  -0.147   0.8838    
## tsWeb:tsDelay -0.005799   0.034988  -0.166   0.8690    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.323 on 52 degrees of freedom
## Multiple R-squared:  0.4935, Adjusted R-squared:  0.4351 
## F-statistic: 8.446 on 6 and 52 DF,  p-value: 2.066e-06
  # 係数を取得して
  agCoef <- coef(oModel)
  
  # 月を行にした係数
  dfEffect <- tibble( nMonth =  1:12) %>%
    mutate(
      gEffect_TV  = agCoef["tsTV"]  + agCoef["tsTV:tsDelay"] * c(9:11, 0:8)[nMonth],   
      gEffect_Web = agCoef["tsWeb"] + agCoef["tsWeb:tsDelay"] * c(9:11, 0:8)[nMonth], 
      gTotal_TV = gEffect_TV / (1-agCoef['L(tsSales)']),
      gTotal_Web = gEffect_Web / (1-agCoef['L(tsSales)'])
    )
  
  # 月xラグを行にした係数
  dfAR <- tibble( nLag = 0:12 ) %>%
    mutate( gAR = agCoef['L(tsSales)']^nLag ) 
  dfEffect2 <- crossing(nMonth = 1:12, nLag = 0:12) %>%
    left_join(dfEffect, by = "nMonth") %>%
    left_join(dfAR, by = "nLag") %>%
    mutate(
      gEffect2_TV = gEffect_TV * gAR, 
      gEffect2_Web = gEffect_Web * gAR
    ) 
  
  # 月xラグのチャート
  dfPlot <- dfEffect2 %>%
    dplyr::select(nMonth, nLag, gEffect2_TV, gEffect2_Web) %>%
    gather(sVar, gValue, starts_with("gEffect2")) %>%
    arrange(nMonth) %>%
    mutate(
      fMonth = factor(nMonth, levels=1:12, labels=paste0(1:12, "月")), 
      fVar = factor(sVar, levels = c("gEffect2_TV", "gEffect2_Web"), labels = c("TV", "Web"))
    )
  g <- ggplot(data = dfPlot, aes(x = nLag, y = gValue, group = fMonth, color=fMonth))
  g <- g + geom_line()
  g <- g + facet_grid(fVar ~ .)
  g <- g + scale_x_continuous(breaks = seq(0, 12))
  g <- g + labs(color = NULL, x = "ラグ", y = "効果")
  g1 <- g + theme_bw()

  # 月のチャート
  dfPlot <- dfEffect %>%
    dplyr::select(nMonth, gTotal_TV, gTotal_Web) %>%
    gather(sVar, gValue, -nMonth) %>%
    mutate(
      fVar = factor(sVar, levels = c("gTotal_TV", "gTotal_Web"), labels = c("TV", "Web"))
    )
  g <- ggplot(data = dfPlot, aes(x = nMonth, y = gValue, fill=fVar))
  g <- g + geom_col()
  g <- g + scale_x_continuous(breaks = seq(1, 12))
  g <- g + labs(x = "月", y = "効果の合計")
  g <- g + guides(fill = FALSE)
  g <- g + facet_grid(fVar ~ .)
  g2 <- g + theme_bw()
  g <- g1 + g2
  g 

  ggsave(paste0(WORKDIR, "/chap8_code1.png"), width = 24, height = 15, units = "cm")

Code 2. 状態空間モデルで部分調整モデルの係数を時間変動させる

  agT <- array(0, dim = c(5, 5, nrow(dfSalesAds)))
  agT[1, 1, ] <- 0.8 # lambda 
  agT[2, 2, ] <- 0.8 # lambda
  agT[3, 3, ] <- 0.8 # lambda
  agT[2, 4, ] <- c(dfSalesAds$gTV[-1], 0)
  agT[3, 5, ] <- c(dfSalesAds$gWeb[-1], 0)
  agT[4, 4, ] <- 1
  agT[5, 5, ] <- 1

  oModel <- SSModel(
    dfSalesAds$gSales ~  -1 +
      SSMcustom(
        Z = matrix(c(1, 1, 1, 0, 0), nrow=1), 
        T = agT,
        R = matrix(c(
          1,   0,  0, 
          0,   0,  0,
          0,   0,  0, 
          0,   1,  0, 
          0,   0,  1
        ), byrow = 1, nrow = 5), 
        Q = diag(NA, nrow = 3) , 
        a1 = c(mean(dfSalesAds$gSales), mean(dfSalesAds$gTV), mean(dfSalesAds$gWeb), 0.5, 0.5),
        # P1 = diag(c(10000, 10000, 10000, 10000, 10000)),
        # P1inf = diag(c(0,0,0,0,0)),
        state_names = c("mu", "stock1", "stock2", "beta1", "beta2")
      ), 
    H = matrix(NA) 
  )
  # # 推定できることを確認
  # oFitted <- fitSSM(
  #   oModel,
  #   inits  = log(c(0.01, 0.01, 0.01, 0.01))#,
  #   #method = "BFGS"
  # )
  # print(oFitted$model$T[,,1])
  # oEstimated <- KFS(oFitted$model)
  # print(summary(oEstimated$alphahat))
  # stop()
  
  # 「パラメータを与えるとモデルを返す」関数を定義する
  sub_update <- function(par, model) {
    # par: 順に{H, sigma^2_epsilon, sigma^2_v1, sigma^2_v2, lambda}
    # model: 現在のモデル
    
    gMag <- 0.5 # 分散の最適値をどのくらい細かく探すか
    # print(c(gMag*exp(par[1:4]), artransform(par[5])))
    newmodel <- model
    newmodel$H[1,1, ] <- exp(gMag*par[1])
    newmodel$Q[1,1, ] <- exp(gMag*par[2])
    newmodel$Q[2,2, ] <- exp(gMag*par[3])
    newmodel$Q[3,3, ] <- exp(gMag*par[4])
    newmodel$T[1,1, ] <- artransform(par[5])
    newmodel$T[2,2, ] <- artransform(par[5])
    newmodel$T[3,3, ] <- artransform(par[5])
    return(newmodel)
  }
  # 推定
  oFitted   <- fitSSM(
    oModel,
    updatefn = sub_update,
    inits  = c(log(0.01), log(0.01), log(0.01), log(0.01), log(0.01), 0.8)
  )
  cat("T: \n")
## T:
  print(oFitted$model$T[,,1])
##        mu stock1 stock2 beta1 beta2
## mu     -1      0      0     0   0.0
## stock1  0     -1      0    12   0.0
## stock2  0      0     -1     0   6.4
## beta1   0      0      0     1   0.0
## beta2   0      0      0     0   1.0
  cat("Q: \n")
## Q:
  print(oFitted$model$Q[,,1])
##          [,1]     [,2]     [,3]
## [1,] 26.92171 0.000000 0.000000
## [2,]  0.00000 0.096265 0.000000
## [3,]  0.00000 0.000000 0.117561
  cat("H: \n")
## H:
  print(oFitted$model$H[,,1])
## [1] 32.84269
  oEstimated <- KFS(oFitted$model)
  print(summary(oEstimated$alphahat))
##        mu              stock1             stock2            beta1      
##  Min.   :-21.427   Min.   :-7.39004   Min.   :-6.7321   Min.   :0.500  
##  1st Qu.:-11.268   1st Qu.: 0.09836   1st Qu.:-0.6677   1st Qu.:2.724  
##  Median :  1.288   Median :12.10449   Median : 2.9516   Median :2.888  
##  Mean   :  0.568   Mean   :11.38453   Mean   : 3.0625   Mean   :2.732  
##  3rd Qu.: 12.459   3rd Qu.:21.46879   3rd Qu.: 6.2365   3rd Qu.:3.031  
##  Max.   : 22.091   Max.   :30.35146   Max.   :15.8279   Max.   :3.353  
##      beta2      
##  Min.   :0.500  
##  1st Qu.:2.147  
##  Median :2.317  
##  Mean   :2.273  
##  3rd Qu.:2.795  
##  Max.   :2.848
  oEstimated <- KFS(oFitted$model)
  print(oEstimated$alphahat)
## Time Series:
## Start = 1 
## End = 60 
## Frequency = 1 
##              mu      stock1     stock2     beta1     beta2
##  1  16.36195055  8.45833333  2.8533333 0.5000000 0.5000000
##  2   1.15350631 -2.45833333  0.3466667 0.8525419 0.6359286
##  3  -3.34713118 11.58053195  1.6883048 1.2890128 0.8025101
##  4  10.24949747 -2.81524479 -1.6883048 1.5576525 0.9690917
##  5  -9.41323102 12.62845575  4.9832166 1.8074535 1.1232572
##  6  12.69675734 -0.69926233 -2.8490279 1.9797638 1.2501797
##  7 -11.86700787 20.69487630  4.8493155 2.1221076 1.3713049
##  8  12.40560755  0.10177824  0.6359042 2.2455777 1.4830224
##  9 -11.57078595 21.68032560  3.5165583 2.3400923 1.5845324
## 10  13.17719152 -0.38548548 -1.9320259 2.3823357 1.6790277
## 11 -11.30776981 18.25300356  6.2974980 2.3744449 1.7522983
## 12   9.77697704  9.29055783 -1.7415224 2.4300492 1.8429489
## 13 -11.51710453 16.46796407  6.9017794 2.5516093 1.9548760
## 14  14.90108763 -2.94443492 -0.2552009 2.6090380 2.0165610
## 15 -15.60500044 25.12125834  5.2966035 2.6878614 2.0859306
## 16  14.89372435 -4.96229789  2.4213398 2.7857598 2.1667923
## 17 -11.77879836 21.39828051  3.6456787 2.8179430 2.2095679
## 18  11.82242624 -1.95447384  3.6458955 2.8490498 2.2517149
## 19 -11.25504116 18.19405776  6.7119929 2.8685923 2.2824646
## 20  12.64607390 -1.55622213  2.4178657 2.8592859 2.2889171
## 21 -11.78399825 18.99786592  3.7622105 2.8311758 2.2852055
## 22  14.28254622 -7.39004518  6.5212141 2.7664357 2.2323962
## 23 -16.27838255 30.35146134  4.4175271 2.7609293 2.2222922
## 24  15.31777542 -6.60746973  5.1383291 2.7849628 2.2302256
## 25 -12.90135556 23.87423907  2.8904829 2.7554253 2.2001720
## 26  12.62046769  3.68001424  1.9498954 2.7359317 2.1728168
## 27 -14.93702429 23.13211622  8.2623436 2.7976153 2.1930063
## 28  18.46556451  0.08809118 -2.7798279 2.7545768 2.1746749
## 29 -20.90549593 27.73313437 15.8278772 2.7996562 2.2202713
## 30  20.53956175  2.50315305 -2.2842225 2.8588673 2.2756151
## 31 -21.42740370 28.94438762 10.2488755 2.9530000 2.3445285
## 32  22.09148050  2.94801260 -6.7320827 3.0214874 2.4090921
## 33 -21.26701174 28.47545642 10.1048117 3.0593147 2.4686153
## 34  20.84103812 -0.94162401 -3.1926887 3.1108505 2.5333469
## 35 -19.21665498 25.82842841  9.5260560 3.1159194 2.5803452
## 36  19.25765613 -2.77062457 -3.8492966 3.1199034 2.6269496
## 37 -16.29664551 23.36198710  7.0016361 3.0540079 2.6580379
## 38  14.67044512 -1.37313010 -0.8881488 3.0299795 2.7054591
## 39 -10.40622386 17.73501919  4.9463375 2.9236132 2.7249490
## 40   8.21080905 15.88653279 -4.9463375 2.9075246 2.7444389
## 41 -10.25635561 19.00376241  8.7885520 2.9792081 2.7764342
## 42  10.82580959  3.63821880  1.7618979 3.0354162 2.7989801
## 43  -8.04682026 18.21677819  4.6757563 3.0200785 2.7936150
## 44   7.58911481  2.92377121  5.3812577 3.0161972 2.7954452
## 45  -7.00620052 25.12686243 -3.4244460 2.9929314 2.7954936
## 46   6.91710269  0.61234763  9.2949826 2.9724055 2.7963591
## 47  -7.47481467 25.84206147  0.2126383 2.9696283 2.8055049
## 48   7.36286216 -0.60022094  9.6066289 2.9702537 2.8163618
## 49  -6.19823945 25.25332695 -0.5942712 2.9363148 2.8109446
## 50   3.52075784  7.33976702  6.2161604 3.0086469 2.8289113
## 51  -0.07265751 19.73805520 -2.5385757 2.9700135 2.8273038
## 52  -3.34258655 14.71410158  5.0831492 3.0730395 2.8391186
## 53   6.46820812  5.56795936  1.7307355 3.1023013 2.8181760
## 54  -6.81607701 21.73229196  1.9328934 3.1425093 2.7992083
## 55   6.81595520  5.92178946  2.8257607 3.1827210 2.7802414
## 56  -5.20080438 22.40442789  3.0127463 3.1715322 2.7464633
## 57   1.42204894  7.09082134 11.8181553 3.2860034 2.8017903
## 58   2.02237736 18.21140496 -4.5335005 3.3056386 2.8180107
## 59  -3.07458142  6.25032091 13.2693337 3.3531157 2.8484748
## 60   5.32501260 13.53306144 -4.4390618 3.3531157 2.8484748
  dfResult <- tibble(
    dMonth = dfSalesAds$dMonth,
    gBeta1 = as.vector(oEstimated$alphahat[,"beta1"]),
    gBeta2 = as.vector(oEstimated$alphahat[,"beta2"]),
    gEffect1 = gBeta1 / (1-oEstimated$model$T[1,1,1]),
    gEffect2 = gBeta2 / (1-oEstimated$model$T[1,1,1])
  )
  print(summary(dfResult))
##      dMonth               gBeta1          gBeta2         gEffect1    
##  Min.   :2014-01-01   Min.   :0.500   Min.   :0.500   Min.   :0.250  
##  1st Qu.:2015-03-24   1st Qu.:2.724   1st Qu.:2.147   1st Qu.:1.362  
##  Median :2016-06-16   Median :2.888   Median :2.317   Median :1.444  
##  Mean   :2016-06-16   Mean   :2.732   Mean   :2.273   Mean   :1.366  
##  3rd Qu.:2017-09-08   3rd Qu.:3.031   3rd Qu.:2.795   3rd Qu.:1.516  
##  Max.   :2018-12-01   Max.   :3.353   Max.   :2.848   Max.   :1.677  
##     gEffect2    
##  Min.   :0.250  
##  1st Qu.:1.073  
##  Median :1.158  
##  Mean   :1.137  
##  3rd Qu.:1.398  
##  Max.   :1.424
  dfPlot <- dfResult %>%
    dplyr::select(dMonth, gEffect1, gEffect2) %>%
    gather(sVar, gValue, -dMonth) %>%
    mutate(fVar = factor(sVar, levels=c("gEffect1", "gEffect2"), labels = c("TV", "Web")))
  g <- ggplot(data=dfPlot, aes(x = dMonth, y = gValue, fill = fVar))
  g <- g + geom_col()
  g <- g + facet_grid(fVar ~ .)
  g <- g + labs(x = "年月", y = "効果の合計")
  g <- g + guides(fill = F)
  g <- g + theme_bw()
  print(g)

  ggsave(paste0(WORKDIR, "/chap8_code2.png"), width = 24, height = 15, units = "cm")

以上