#----------------------------
# 関数の読み込み（plot1_1.R）
#----------------------------

source("plot1_1.R")　#図1-4のコードも参照

#------------------
# データの読み込み
#------------------

dat1_1 <- read.csv("approve_economy.csv", fileEncoding = "SJIS") %>%
  mutate(
    sociotropic_gd = q5c1 + q5c2,
    sociotropic_bd = q5c4 + q5c5,
    liv_gd         = q3c1 + q3c2,
    liv_bd         = q3c4 + q3c5
  )

df_keiki <- dat1_1 %>%
  select(yearmon, year_graph, sociotropic_gd, sociotropic_bd)

df_kurashi <- dat1_1 %>%
  select(yearmon, year_graph, liv_gd, liv_bd)

#------------------
#  プロットの作成
#------------------

a1_1 <- plot1_1(df_keiki, "sociotropic_bd", "sociotropic_gd", "景気への評価(-)", "景気への評価(+)")
a1_2 <- plot1_1(df_kurashi, "liv_bd", "liv_gd", "暮らし向きへの評価(-)", "暮らし向きへの評価(+)")


grid.arrange(a1_1[[1]], a1_1[[2]], a1_1[[3]], ncol = 3, widths = c(4/9, 1/9, 4/9))
grid.arrange(a1_2[[1]], a1_2[[2]], a1_2[[3]], ncol = 3, widths = c(4/9, 1/9, 4/9))

#---------------------------
# 関数の読み込み(plot1_3.R)
#---------------------------

source("plot1_3.R")

#------------------
# データの読み込み
#------------------

dat1_3 <- read.csv("df_ec_keiki.csv") %>%
  mutate(否定 = 否定*100,
         肯定 = 肯定*100)
  
#----------------
# プロットの作成
#----------------

a1_3 <- plot1_3(dat1_3, "景気への評価(-)", "景気への評価(+)")

grid.arrange(a1_3[[1]], a1_3[[2]], a1_3[[3]], ncol = 3, widths = c(4/9, 1/9, 4/9))
#景気評価に関する図は出力を省略。

#------------------
# データの読み込み
#------------------

dat1_4 <- read_csv("approve_economy.csv", locale = locale(encoding = "SJIS")) %>%
  dplyr::select(yearmon, year_graph, ruling, opposite, indep)
#（上述の"approve_economy.csv"データに依拠、データの閲覧・入手については筆者にご相談ください）

#---------------------------------------
# プロットの作成(plot1_1.Rの内部も参照)
#---------------------------------------

a.mid<-ggplot(dat1_4,aes(x=1,y=year_graph))+
  geom_text(aes(label=year_graph))+
  geom_segment(aes(x=0.94,xend=0.96,yend=year_graph))+
  geom_segment(aes(x=1.04,xend=1.065,yend=year_graph))+
  ggtitle("")+
  ylab(NULL)+
  theme(axis.title=element_blank(),
        panel.grid=element_blank(),
        axis.text.y=element_blank(),
        axis.ticks.y=element_blank(),
        panel.background=element_blank(),
        axis.text.x=element_text(color=NA),
        axis.ticks.x=element_line(color=NA),
        plot.margin = unit(c(1,-1,1,-1), "mm")) +
  scale_x_reverse()


a1 <- ggplot(data = dat1_4, aes(x = as.Date(yearmon,"%Y-%m-%d"), y = ruling)) +
  geom_bar(stat = "identity") +  ggtitle("与党支持率") +
  theme(axis.title.x = element_blank(), 
        axis.title.y = element_blank(), 
        axis.text.y = element_blank(), 
        axis.ticks.y = element_blank(), 
        plot.margin = unit(c(1,-1,1,0), "mm")) +
  coord_flip()+
  ylim(0, 75)+
  xlab("年")+
  ylab("与党支持率")+
  theme(text = element_text(family = "jpfont"))+
  theme_bw()


a2 <- ggplot(data = dat1_4, aes(x = as.Date(yearmon,"%Y-%m-%d"), y = opposite)) +
  geom_bar(stat = "identity") + ggtitle("野党支持率") +
  theme(axis.title.x = element_blank(), 
        axis.title.y = element_blank(), 
        axis.text.y = element_blank(), 
        axis.ticks.y = element_blank(), 
        plot.margin = unit(c(1,-1,1,0), "mm")) +
  coord_flip()+
  ylim(0, 75)+
  xlab("年")+
  ylab("野党支持率")+
  theme(text = element_text(family = "jpfont"))+
    theme_bw()


a3 <- ggplot(data = dat1_4, aes(x = as.Date(yearmon,"%Y-%m-%d"), y = indep)) +
  geom_bar(stat = "identity") + ggtitle("無党派層割合") +
  theme(axis.title.x = element_blank(), 
        axis.title.y = element_blank(), 
        axis.text.y = element_blank(), 
        axis.ticks.y = element_blank(), 
        plot.margin = unit(c(1,-1,1,0), "mm")) +
  coord_flip()+
  xlab("年")+
  ylab("無党派割合")+
  ylim(0, 75)+
  theme(text = element_text(family = "jpfont"))+
    theme_bw()


aa1 <- ggplot_gtable(ggplot_build(a1))
aa2 <- ggplot_gtable(ggplot_build(a2))
aa3 <- ggplot_gtable(ggplot_build(a3))

aa.mid <- ggplot_gtable(ggplot_build(a.mid))

grid.arrange(aa1,aa.mid,aa2,aa.mid,aa3,ncol=5,widths=c(4/14,1/14,4/14,1/14,4/14))

#------------------
# データの読み込み
#------------------

dat1_5 <- read_csv("psupport_move.csv")

#------------
# 要件の設定
#------------

dat1_5 <- dat1_5 %>%
  mutate(Condition = factor(Condition, levels = c("残留: 与党→与党", "残留: 野党→野党", 
                                                  "離脱: 与党→野党", "離脱: 野党→与党")))

#----------------
# プロットの作成
#----------------

ggplot(dat1_5, aes(x = 年, y = Ratio * 100, fill = Condition)) +
  geom_bar(stat = "identity", position = position_dodge(), color = "black") +
  geom_text(aes(label = sprintf("%.1f", Ratio * 100)), 
            position = position_dodge(width = 0.9), 
            vjust = -0.5, size = 4, color = "black") +
  facet_wrap(~Condition, ncol = 2) +
  scale_fill_grey(start = 0.5, end = 0.8) + 
  labs(title = "", x = "年", y = "割合(%) ") +
  theme_bw(base_size = 14) +
  theme(
    legend.position = "none", 
    panel.grid.major = element_line(size = 0.5),
    panel.grid.minor = element_blank(),
    strip.text = element_text(face = "bold", size = 12) 
  ) +
  geom_label(aes(label = sprintf("%.1f", Ratio * 100)), 
             position = position_dodge(width = 0.9), 
             vjust = -0.5, size = 4, fill = "white", color = "black", label.size = 0.3, 
             label.padding = unit(0.2, "lines"))

#------------------
# データの読み込み
#------------------

jes2019 <- read_csv("2019jes.csv", locale = locale(encoding = "SJIS"))

#-------------------------------
# 「与党―与党」のダミー変数作成
#-------------------------------

jes2019 <- jes2019 %>%
  mutate(
    rul_rul_dum2019 = case_when(
      (Q2 %in% c("自民党", "公明党")) & (aQ7 %in% c("自民党", "公明党")) ~ 1,
      TRUE ~ 0
    )
  )

#-------------------------------
# 「与党―野党」のダミー変数作成
#-------------------------------

jes2019 <- jes2019 %>%
  mutate(
    rul_op_dum2019 = case_when(
      (aQ7 %in% c("自民党", "公明党")) &
      (Q2 %in% c("立憲民主", "社民党", "共産党", "次世代の", "みんなの",
                 "希望の党", "新党改革", "日本未来の党", "日本維新", "その他")) ~ 1,
      TRUE ~ 0
    )
  )

library(tidyverse)

#------------------
# データの読み込み
#------------------

dat1_6 <- read_csv("df_temp.csv", locale = locale(encoding = "SJIS"))

#---- 与党派：社会志向(+)
dat1_6_1 <- dat1_6 %>%
  select(yearmon, sociotropic_gd_rul, sociotropic_bd_rul) %>%
  mutate(政党支持 = "与党派")

#---- 野党派：社会志向(+)
dat1_6_2 <- dat1_6 %>%
  select(yearmon, sociotropic_gd_op, sociotropic_bd_op) %>%
  mutate(政党支持 = "野党派")

#---- 無党派：社会志向(+)
dat1_6_3 <- dat1_6 %>%
  select(yearmon, sociotropic_gd_indep, sociotropic_bd_indep) %>%
  mutate(政党支持 = "無党派")

colnames(dat1_6_1) <- colnames(dat1_6_2) <- colnames(dat1_6_3) <- c("年", "肯定", "否定", "政党支持")
dat1_6 <- bind_rows(dat1_6_1, dat1_6_2, dat1_6_3) %>%
  mutate(across(c(肯定, 否定), ~ . * 100))

dat1_6 <- dat1_6 %>%
  pivot_longer(cols = c(肯定, 否定), names_to = "景気評価", values_to = "値") %>%
  mutate(政党支持 = factor(政党支持, levels = c("与党派", "野党派", "無党派")))

#------------------
# プロットの作成
#------------------

ggplot(data = dat1_6, aes(x = as.Date(年,"%Y-%m-%d"), y = 値)) +
  geom_line(aes(color = 景気評価), size = 1) +
  scale_color_manual(values = c("grey", "black")) +
  xlab("年")+
  facet_wrap(~政党支持,ncol=1)+
  ylab("割合(%)")+ 
  theme_bw()+
theme(text = element_text(family = "jpfont"))

#----------------------------
#　関数の読み込み(plot1_7.R)
#----------------------------

source("plot1_7.R")

#------------------------------------------------
# データの読み込みとプロットの作成: 社会志向評価
#------------------------------------------------

dat1_7socio <- read.csv("df_ec_partisan.csv")
a1_7socio <- plot1_7(dat1_7socio, "景気評価", "景気")


grid.arrange(a1_7socio$g1,  a1_7socio$g2, 
              a1_7socio$g3, 
             ncol = 3, widths = c(1/3, 1/3, 1/3))

#------------------
# モジュールの導入
#------------------


import statsmodels.api as sm
import pandas

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

import pandas_datareader as pdr
from statsmodels.regression.rolling import RollingOLS
import matplotlib.pyplot as plt
import seaborn
seaborn.set_style('darkgrid')
pd.plotting.register_matplotlib_converters()
%matplotlib inline


#------------------
# データの読み込み
#------------------

df = pd.read_csv('df_fecm.csv',
                  index_col='yearmon', parse_dates=True,encoding='cp932')
print(df)

#------------------
# 内生変数の設定
#------------------

endog = df.approve
endog

#---------------------------------------
# 外生変数の設定：社会志向の経済評価(+)
#---------------------------------------

exog = sm.add_constant(df['ruling'])

#----------------------
# 推定と推定結果の確認
#----------------------

rols_ruling = RollingOLS(endog, exog, window=60)
rres_ruling = rols_ruling.fit()
params = rres_ruling.params
print(params.head())
print(params.tail())

#------------------
# プロットの作成
#------------------

fig_ruling= rres_dis_liv.plot_recursive_coefficient(variables=['ruling'], figsize=(14,6))

#---------------------------
# 関数の読み込み(plot2_2.R)
#---------------------------
source("plot2_2.R")

#------------------
# データの読み込み
#------------------
dat2_2 <- read.csv("approve_economy.csv", fileEncoding = "SJIS")


# 暮らし向き
cols1 <- c("q3c1", "q3c2", "q3c3", "q3c4", "q3c5", "q3c6")
labels1 <- c("大変楽になった", "やや楽になった", "どちらでもない", "やや苦しくなった", "大変苦しくなった", "わからない")
b2_2_1 <- plot2_1(dat2_2, cols1, labels1, "暮らし向き")

# 物価
cols2 <- c("q4c3", "q4c1", "q4c2", "q4c4")
labels2 <- c("今より下がると思う", "落ち着いてきたと思う", "上がると思う", "わからない")
b2_2_2 <- plot2_1(dat2_2, cols2, labels2, "物価")

# 景気
cols3 <- c("q5c1", "q5c2", "q5c3", "q5c4", "q5c5", "q5c6")
labels3 <- c("確かに良くなってきたと思う", "やや良くなってきたと思う", "どちらでもない", "やや悪くなってきたと思う", "確かに悪くなってきたと思う", "わからない")
b2_2_3 <- plot2_1(dat2_2, cols3, labels3, "景気")

# これから先の生活
cols4 <- c("q6c1", "q6c2", "q6c3", "q6c4")
labels4 <- c("良くなっていく", "変わりない", "悪くなっていく", "どちらともいえない")
b2_2_4 <- plot2_1(dat2_2, cols4, labels4, "これから先の生活")

#------------------
# プロットの作成
#------------------

gridExtra::grid.arrange(b2_2_1, b2_2_2, b2_2_3, b2_2_4, ncol = 2, nrow = 2)

#------------------
# データの読み込み
#------------------

dat2_2 <- read.csv("approve_economy.csv", fileEncoding = "SJIS")


#--------------------------------- 
# FIの計算：社会志向の経済評価(+)
#--------------------------------- 

fdGPH(na.omit(sociotropic_gd))

#--------------------------------- 
# FIの計算：社会志向の経済評価(-)
#--------------------------------- 

fdGPH(na.omit(sociotropic_bd))

#--------------------------------- 
# FIの計算：個人志向の経済評価(+)
#--------------------------------- 

fdGPH(na.omit(liv_gd))

#--------------------------------- 
# FIの計算：個人志向の経済評価(-)
#--------------------------------- 

fdGPH(na.omit(liv_bd))

#---------------- 
# 関数の読み込み
#----------------
source("adf_kpss_test.R")


#------------------------
# データの読み込みと設定
#------------------------ 
dat2_2 <- read.csv("approve_economy.csv", fileEncoding = "SJIS")

vars <- dat2_2 %>%
  dplyr::select(sociotropic_gd, sociotropic_bd, liv_gd, liv_bd) %>%
  as.list() %>%
  setNames(c("社会志向(+)", "社会志向(-)", "個人志向(+)", "個人志向(-)"))


#---------------- 
# 分析結果の表化
#----------------
df_results <- map_dfr(
  names(vars),
  ~ adf_kpss_test(vars[[.x]]) %>%
    mutate(変数名 = .x)   
) %>%
  select(変数名, everything())  

#----------- 
# 表の出力
#----------- 
tb2_2 <- df_results %>%
  kbl(
    format = "html",
    digits = 4,                     
  ) %>%
  kable_styling(full_width = FALSE)

#---------------------------- 
# 関数の読み込み(plot2_3.R)
#---------------------------- 
source(plot2_3.R)

#------------------
# データの読み込み
#------------------
dat2_3 <- read.csv("approve_economy.csv", fileEncoding = "SJIS") %>%
  mutate(
    sociotropic_gd = q5c1 + q5c2,
    sociotropic_bd = q5c4 + q5c5,
    liv_gd         = q3c1 + q3c2,
    liv_bd         = q3c4 + q3c5
  )

#---------------- 
# イベントの設定
#----------------
events <- tibble(
  date = as.Date(c("1973-10-01", "1989-04-01", "1990-01-01", "1997-07-01", 
                   "2008-09-01", "2011-03-01", "2012-12-01", "2020-04-01")),
  label = c("第1次石油ショック(1973)", "消費税導入(1989)", "バブル経済崩壊(1990)",
            "アジア通貨危機(1997)", "グローバル経済危機(2008)", "東日本大震災(2011)",
            "アベノミクス導入(2012)", "緊急事態宣言(2020)")
)

#-------------- 
# データの整形
#--------------
dat2_3keiki <- dat2_3 %>%
  select(yearmon, sociotropic_gd, sociotropic_bd) %>%
  pivot_longer(cols = -yearmon, names_to = "経済評価", values_to = "value") %>%
  mutate(経済評価 = recode(経済評価, sociotropic_gd = "肯定", sociotropic_bd = "否定"))

dat2_3kurashi <- dat2_3 %>%
  select(yearmon, liv_gd, liv_bd) %>%
  pivot_longer(cols = -yearmon, names_to = "経済評価", values_to = "value") %>%
  mutate(経済評価 = recode(経済評価, liv_gd = "肯定", liv_bd = "否定"))

#-----------------
# プロットの作成
#-----------------
b1 <- plot2_3(dat2_3keiki , "A:社会志向の経済評価の推移と政治経済的出来事", 45)
b2 <- plot2_3(dat2_3kurashi, "B:個人志向の経済評価の推移と政治経済的出来事", 28)

windows(30, 20)
grid.arrange(b1, b2, nrow = 2)

#---------------------------- 
# 関数の読み込み(plot2_4.R)
#---------------------------- 
source("plot2_4.R")

#------------------
# データの読み込み
#------------------
dat2_4 <- read.csv("approve_economy.csv", fileEncoding = "SJIS") %>%
  mutate(
    sociotropic_gd = q5c1 + q5c2,
    sociotropic_bd = q5c4 + q5c5,
    liv_gd         = q3c1 + q3c2,
    liv_bd         = q3c4 + q3c5
  )

#---------------------------- 
# 景気循環のセグメントを定義
#----------------------------
business_cycles <- data.frame(
  start = as.Date(c("1963-12-01", "1965-10-01", "1971-12-01", "1975-03-01",
                    "1977-10-01", "1983-02-01", "1986-11-01", "1993-10-01",
                    "1999-12-01", "2002-01-01", "2009-03-01", "2012-11-01")),
  end = as.Date(c("1964-10-01", "1970-07-01", "1973-11-01", "1977-01-01",
                  "1980-02-01", "1985-06-01", "1991-02-01", "1997-05-01",
                  "2000-11-01", "2008-02-01", "2012-03-01", "2020-05-01"))
)



#----------------
# プロットの作成
#----------------
b1 <- plot2_4(dat2_4, "sociotropic_gd", "sociotropic_bd", 
              "A: 社会志向の経済評価の推移と景気循環", 48, "社会志向の経済評価")

b2 <- plot2_4(dat2_4, "liv_gd", "liv_bd", 
              "B: 個人志向の経済評価の推移と景気循環", 18, "個人志向の経済評価")

windows(width = 240, height = 180)
grid.arrange(b1, b2, nrow = 2)

#------------------ 
# データの読み込み
#------------------
dat2_3 <- read_csv("approve_economy.csv", locale = locale(encoding = "SJIS")) %>%
  mutate(
    与党派_社会志向_否定 = sociotropic_bd_rul,
    野党派_社会志向_否定 = sociotropic_bd_op,
    無党派_社会志向_否定 = sociotropic_bd_indep,
    与野党派間差異 = sociotropic_bd_rul - sociotropic_bd_op,
    与無党派間差異 = sociotropic_bd_rul - sociotropic_bd_indep,
    野無党派間差異 = sociotropic_bd_op - sociotropic_bd_indep
  ) %>%
  select(与党派_社会志向_否定, 野党派_社会志向_否定, 無党派_社会志向_否定, 
         与野党派間差異, 与無党派間差異, 野無党派間差異)

#----------------
# 記述統計の抽出
#----------------
stats <- dat2_3 %>%
  describe() %>%
  select(mean, sd, min, max, skew, kurtosis) %>%
  rename(平均 = mean, 標準偏差 = sd, 最小値 = min, 最大値 = max, 歪度 = skew, 尖度 = kurtosis)

#---------- 
# 表の作成
#---------- 
stats %>%
  knitr::kable(
    format = "html",
    booktabs = TRUE, 
    caption = "党派別の社会志向の肯定的経済評価・党派性差異の記述統計",
    digits = 3
  ) %>%
  add_header_above(header = c(" " = 1, "統計量" = 6)) %>%
  kable_styling(bootstrap_options = "striped", full_width = FALSE)

#-----------------
# データの読み込み
#-----------------
dat2_4 <- read_csv("approve_economy.csv", locale = locale(encoding = "SJIS"))

#-------------------
# データリストの作成
#-------------------
dat2_4 <- list(
  dat$sociotropic_gd_rul, dat$sociotropic_gd_op, dat$sociotropic_gd_indep,
  dat$sociotropic_bd_rul, dat$sociotropic_bd_op, dat$sociotropic_bd_indep,
  dat$egotropic_gd_rul,   dat$egotropic_gd_op,   dat$egotropic_gd_indep,
  dat$egotropic_bd_rul,   dat$egotropic_bd_op,   dat$egotropic_bd_indep
)

labels <- c(
  "与党派・社会志向(+)", "野党派・社会志向(+)", "無党派・社会志向(+)",
  "与党派・社会志向(-)", "野党派・社会志向(-)", "無党派・社会志向(-)",
  "与党派・個人志向(+)", "野党派・個人志向(+)", "無党派・個人志向(+)",
  "与党派・個人志向(-)", "野党派・個人志向(-)", "無党派・個人志向(-)"
)

#----------
# FIの計算
#----------
results2_4 <- lapply(data_list2_4, function(data) {
  res <- fdGPH(na.omit(data))
  c(round(res$d, 3), round(res$sd.reg, 3))
})


#------------
# 結果の表化
#------------
results_df2_4 <- results %>%
  do.call(rbind, .) %>%
  as.data.frame() %>%
  setNames(c("d", "sd.reg")) %>%
  mutate(Variable = labels) %>%
  select(Variable, everything())

#-----------------
# データの読み込み
#-----------------
data_list2_4 <- read_csv("approve_economy.csv", locale = locale(encoding = "SJIS")) %>%
  select(
    sociotropic_gd_rul, sociotropic_gd_op, sociotropic_gd_indep,
    sociotropic_bd_rul, sociotropic_bd_op, sociotropic_bd_indep,
    egotropic_gd_rul,   egotropic_gd_op,   egotropic_gd_indep,
    egotropic_bd_rul,   egotropic_bd_op,   egotropic_bd_indep
  ) %>%
  as.list()

#----------------------
#各変数に対応するラベル
#----------------------
labels <- c(
  "与党派・社会志向(+)", "野党派・社会志向(+)", "無党派・社会志向(+)",
  "与党派・社会志向(-)", "野党派・社会志向(-)", "無党派・社会志向(-)",
  "与党派・個人志向(+)", "野党派・個人志向(+)", "無党派・個人志向(+)",
  "与党派・個人志向(-)", "野党派・個人志向(-)", "無党派・個人志向(-)"
)

#----------------
# 分析結果の表化
#----------------
df_results2_4 <- map2_dfr(data_list2_4, labels, ~{
  series <- na.omit(.x)       
  
  # ADF 検定
  adf_res <- adf.test(series)
  adf_stat <- round(adf_res$statistic, 3)
  adf_pval <- round(adf_res$p.value, 4)
  
  # KPSS 検定
  kpss_res <- kpss.test(series, null = "Level")
  kpss_stat <- round(kpss_res$statistic, 3)
  kpss_pval <- round(kpss_res$p.value, 4)
  
  tibble(
    変数名  = .y,         
    ADF_Stat  = adf_stat,
    ADF_pval  = adf_pval,
    KPSS_Stat = kpss_stat,
    KPSS_pval = kpss_pval
  )
})

#----------
# 表の出力
#----------
tb2_4 <- df_results2_4 %>%
  kbl(
    format = "html",
    digits = 3
  ) %>%
  kable_styling(full_width = FALSE)

#------------------ 
# 関数の読み込み
#------------------
source("plot2_5df.R")
source("plot2_5_1.R")
source("plot2_5_2.R")


#------------------ 
# データの読み込み
#------------------
dat2_5 <- read.csv("approve_economy.csv", fileEncoding = "SJIS")
attach(dat2_5)
year <- read.csv("year.csv")
attach(year)


#------------------ 
# ハイライトの設定
#------------------
admin <- data.frame(
  start = as.Date( '2009-09-01'), 
  end   = as.Date('2012-11-01'))


#------------------ 
# プロットの作成
#------------------
#---- 社会志向 (+)：原系列
df_parec1 <- plot2_5df(year$yearmon, cbind(sociotropic_gd_rul, sociotropic_gd_op), c("年", "与党派", "野党派"))
b1 <- plot2_5_1(df_parec1, "A:党派別の社会志向の経済評価(+)","年", "value", "党派別経済効果")

#---- 社会志向 (+): 差異
socposdiff = sociotropic_gd_rul-sociotropic_gd_op
df_parec2 <- create_dataframe(year$yearmon, socposdiff, c("年", "党派性差異"))
b2 <- plot2_5_2(df_parec2, "B:党派別の社会志向の党派性差異(+)", "党派性差異", "value")


grid.arrange(b1, b2, ncol = 1)

#------------------ 
# プロットの作成
#------------------
#---- 社会志向 (-)：2系列併記
df_parec3 <- plot2_5df(year$yearmon, cbind(sociotropic_bd_rul, sociotropic_bd_op), c("年", "与党派", "野党派"))
b3 <- plot2_5_1(df_parec3, "A:党派別の社会志向の経済評価(-)","年", "value", "党派別経済効果")

#---- 社会志向 (-): 差異
socnegdiff = sociotropic_bd_rul-sociotropic_bd_op
df_parec4 <- plot2_5df(year$yearmon, socnegdiff, c("年", "党派性差異"))
b4 <- plot2_5_2(df_parec4, "B:党派別の社会志向の党派性差異(-)", "党派性差異", "value")


grid.arrange(b3, b4, ncol = 1)

#------------------------
# データの読み込みと整形
#------------------------

dat3_1 <- read_csv("sentiment_data_includingNEWS.csv", locale = locale(encoding = "SJIS")) %>%
  na.omit() %>%
  group_by(newspaper) %>%
  summarise(across(c(positive, negative), mean, na.rm = TRUE)) %>%
  ungroup() %>%
  pivot_longer(cols = c(positive, negative), names_to = "sentiment", values_to = "value") %>%
  mutate(
    sentiment = recode(sentiment, "positive" = "肯定", "negative" = "否定"),
    newspaper = recode(newspaper, "asahi" = "朝日", "nikkei" = "日経", "yomiuri" = "読売"),
    newspaper = factor(newspaper, levels = c("朝日", "日経", "読売")),
    sentiment = factor(sentiment, levels = c("肯定", "否定"))
  )

#-----------
# ANOVA検定
#-----------
anova_result <- dat3_1 %>% aov(value ~ newspaper * sentiment, data = .)
summary(anova_result)

#----------------
# Tukey HSD 検定
#----------------
tukey_result <- TukeyHSD(anova_result)
print(tukey_result)


#----------------
# プロットの作成
#----------------
windows(10, 8)
dodge <- position_dodge(width = 0.8)

ggplot(dat3_1, aes(x = newspaper, y = value, fill = sentiment, group = sentiment)) +
    geom_bar(stat = "identity", position = dodge, width = 0.7) +  
    geom_label(aes(label = round(value, 2)),
             position = dodge, vjust = -0.2, size = 3, fill = "white", color = "black", label.size = 0.25) +
  xlab("新聞社") +
  ylab("肯定的／否定的情報スコア") +
  geom_signif(comparisons = list(c("朝日", "日経"), c("日経", "読売"), c("朝日", "読売")),
              map_signif_level = TRUE, y_position = c(0.34, 0.34, 0.36),
              annotation = c("***.", "***", "*"), tip_length = 0) +
  geom_signif(
    y_position = c(0.25, 0.28, 0.32),
    xmin = c(0.8, 1.8, 2.8),
    xmax = c(1.2, 2.2, 3.2),
    annotation = c("NS.", "***", "***"),
    tip_length = 0
  ) +
  theme_bw() +
  scale_fill_grey() +
  labs(fill = "情報")

#------------------------------------
# データの読み込みと度数分布表の作成
#------------------------------------
dat3_2 <- read.csv("newspaper_text.csv", fileEncoding = "SJIS") %>% 
  { table(.$newspaper) } %>%           
  as.data.frame() %>%                  
  { cbind(c("朝日", "日経", "読売"), .) } %>% 
  t()                                   

#------------------
# データの読み込み
#------------------
dat3_3 <- read.csv("approve_economy.csv", fileEncoding = "SJIS")
attach(dat3_3)

#-------------------------
# 経済情報スコアの記述統計
#-------------------------
des_ts <- describe(data.frame(positive, negative))

des_stats <- des_ts[, c("n", "mean", "sd", "median", "min", "max")]
colnames(des_stats) <- c("n", "平均", "標準偏差", "中央値", "最小値", "最大値")

#-------------------------
# 経済情報スコアのFI分析
#-------------------------
#肯定
fdGPH(na.omit(positive))

#否定
fdGPH(na.omit(negative))

#------------------
# データの読み込み
#------------------
dat3_2_1 <- read.csv("approve_economy.csv", fileEncoding = "SJIS") %>%
  dplyr::select(yearmon, positive, negative) %>%
  stats::setNames(c("yearmon", "肯定", "否定")) %>%
  tidyr::gather("経済報道", "value", -yearmon)

#-------------------------------
# A: 経済報道と政治経済的出来事
#-------------------------------
c1 <- dat3_2_1 %>%
  ggplot2::ggplot(aes(x = as.Date(yearmon, "%Y-%m-%d"), y = value)) +
  geom_line(aes(linetype = 経済報道)) +
  theme_minimal() +
  xlab("年") +
  ylab("肯定的／否定的経済報道") +
  # 重要イベントの縦線 & 注釈
  annotate("text", x = as.Date("1976-02-01", "%Y-%m-%d"),
           y = 2, label = "第1次石油ショック(1973)", angle = 45, size = 2) +
  geom_vline(xintercept = as.Date("1973-10-01")) +
  annotate("text", x = as.Date("1989-01-01", "%Y-%m-%d"),
           y = 2, label = "消費税導入(1989)", angle = 45, size = 2) +
  geom_vline(xintercept = as.Date("1989-04-01")) +
  annotate("text", x = as.Date("1991-01-01", "%Y-%m-%d"),
           y = 2, label = "バブル経済崩壊(1990)", angle = 45, size = 2) +
  geom_vline(xintercept = as.Date("1990-01-01")) +
  annotate("text", x = as.Date("1997-09-01", "%Y-%m-%d"),
           y = 2, label = "アジア通貨危機(1997)", angle = 45, size = 2) +
  geom_vline(xintercept = as.Date("1997-07-01")) +
  annotate("text", x = as.Date("2008-11-01", "%Y-%m-%d"),
           y = 2, label = "グローバル経済危機(2008)", angle = 45, size = 2) +
  geom_vline(xintercept = as.Date("2008-09-01")) +
  annotate("text", x = as.Date("2011-03-01", "%Y-%m-%d"),
           y = 2, label = "東日本大震災(2011)", angle = 45, size = 2) +
  geom_vline(xintercept = as.Date("2011-03-01")) +
  annotate("text", x = as.Date("2013-03-01", "%Y-%m-%d"),
           y = 2, label = "アベノミクス導入(2012)", angle = 45, size = 2) +
  geom_vline(xintercept = as.Date("2012-12-01")) +
  annotate("text", x = as.Date("2020-04-01", "%Y-%m-%d"), 
           y = c(45), label = c("緊急事態宣言(2020)"), angle=315,size=5)+
  geom_vline(xintercept = as.Date("2020-04-01"))+
  ggtitle("A:経済報道と政治経済的出来事")


#------------------------
# B: 経済報道と景気循環
#------------------------
dat3_2_2 <- read.csv("approve_economy.csv", fileEncoding = "SJIS") %>%
  dplyr::select(yearmon, positive, negative) %>%
  stats::setNames(c("yearmon", "肯定", "否定")) %>%
  tidyr::gather("news", "value", -yearmon)

c2 <- dat3_2_2 %>%
  ggplot2::ggplot(aes(x = as.Date(yearmon, "%Y-%m-%d"), y = value)) +
  geom_line(aes(linetype = news), size = 1) +
  theme_minimal() +
  xlab("年") +
  ylab("肯定的／否定的経済報道") +
  ggtitle("B:経済報道と景気循環") +
  geom_segment(aes(x=as.Date("1963-12-01"),
                   xend=as.Date("1964-10-01"),y=0.035,yend=0.035,size=1.5))+
  geom_segment(aes(x=as.Date("1965-10-01"),
                   xend=as.Date("1970-07-01"),y=0.035,yend=0.035,size=1.5))+
  geom_segment(aes(x=as.Date("1971-12-01"),
                   xend=as.Date("1973-11-01"),y=0.035,yend=0.035,size=1.5))+
  geom_segment(aes(x=as.Date("1975-03-01"),
                   xend=as.Date("1977-01-01"),y=0.035,yend=0.035,size=1.5))+
  geom_segment(aes(x=as.Date("1977-10-01"),
                   xend=as.Date("1980-02-01"),y=0.035,yend=0.035,size=1.5))+
  geom_segment(aes(x=as.Date("1983-02-01"),
                   xend=as.Date("1985-06-01"),y=0.035,yend=0.035,size=1.5))+
  geom_segment(aes(x=as.Date("1986-11-01"),
                   xend=as.Date("1991-02-01"),y=0.035,yend=0.035,size=1.5))+
  geom_segment(aes(x=as.Date("1993-10-01"),
                   xend=as.Date("1997-05-01"),y=0.035,yend=0.035,size=1.5))+
  geom_segment(aes(x=as.Date("1999-12-01"),
                   xend=as.Date("2000-11-01"),y=0.035,yend=0.035,size=1.5))+
  geom_segment(aes(x=as.Date("2002-01-01"),
                   xend=as.Date("2008-02-01"),y=0.035,yend=0.035,size=1.5))+
  geom_segment(aes(x=as.Date("2009-03-01"),
                   xend=as.Date("2012-03-01"),y=0.035,yend=0.035,size=1.5))+
  geom_segment(aes(x=as.Date("2012-11-01"),
                   xend=as.Date("2020-05-01"),y=0.035,yend=0.035,size=1.5))+
  theme(legend.position = "none")


#----------------
# プロットの作成
#----------------
windows(30, 20)
gridExtra::grid.arrange(c1, c2, nrow = 2)

*************************************************************
* 経済評価についての実数和分誤差修正メカニズムの分析（FECM）*
*************************************************************

/*
必要に応じて、以下を実行
 clear
 cd "C:Users\yourpath"
 use " ", clear
 */

***---time series setting
tsset time_id

***----variables setting
drop rcpi_socgd-_est_reg39 
*clear後であれば不要

gen sociotropic_gd=q5c1+q5c2
gen sociotropic_bd=q5c4+q5c5

gen liv_gd=q3c1+q3c2
gen liv_bd=q3c4+q3c5

gen ln_nav = log(nav_high)
gen lnIIP = log(iip)


*************************
* 社会志向の肯定的評価  *
*************************

***---1st step: 残差の予測値の計算----*


reg sociotropic_gd cpi
predict rcpi_socgd, resid

reg sociotropic_gd unemp
predict runemp_socgd, resid

reg sociotropic_gd ln_nav
predict rnav_socgd, resid

reg sociotropic_gd ruling
predict rldp_socgd, resid

reg sociotropic_gd m2
predict rm2_socgd, resid

reg sociotropic_gd iip
predict riip_socgd, resid

reg sociotropic_gd interest
predict rint_socgd, resid

reg sociotropic_gd positive
predict rpositive_socgd, resid

reg sociotropic_gd negative
predict rnegative_socgd, resid


reg sociotropic_gd ruling positive negative  cpi unemp ln_nav lnIIP   business_cycle crisis_dum 
predict rsocgd_full, resid

reg sociotropic_gd  positive negative  cpi unemp ln_nav lnIIP  business_cycle crisis_dum 
predict rsocgd_no, resid


***---2nd step: d値によるフィルタリンリングの実施 ----*

arfima sociotropic_gd,  iterate(2000)
predict d_socgd, fdifference

arfima liv_gd, iterate(2000)
predict d_livgd, fdifference

arfima cpi, iterate(2000)
predict d_cpi, fdifference

arfima unemp, iterate(2000)
predict d_unemp, fdifference

arfima ln_nav, iterate(2000)
predict d_nav, fdifference

arfima ruling, iterate(2000)
predict d_rul, fdifference

arfima iip, iterate(2000)
predict d_iip, fdifference


arfima interest, iterate(2000)
predict d_int, fdifference

arfima positive, iterate(2000)
predict d_positive, fdifference

arfima negative, iterate(2000)
predict d_negative, fdifference

arfima indep, iterate(2000)
predict d_indep, fdifference


arfima rsocgd_full, iterate(2000)
predict d_rsocgd_full, fdifference


arfima rsocgd_no, iterate(2000)
predict d_rsocgd_no, fdifference


***---3rd step: OLSによるFECM推定----*

* 各種のパラメターの設定
 scalar ecmsig =-1.645 
 scalar ecmMCV =-3.867 
 scalar up = 1.96 
 scalar lp =-1.96 
 scalar alpha = .05

*******党派性無モデル********

*----AICの比較から選択される特定化
reg d_socgd  d_positive d_negative l.rpositive_socgd
*Akaike information criteria
estat ic
*\
Akaike's information criterion and Bayesian information criterion

-----------------------------------------------------------------------------
       Model |        Obs  ll(null)  ll(model)      df         AIC        BIC
-------------+---------------------------------------------------------------
           . |        656 -1511.357  -1476.208       4    2960.417   2978.361
-----------------------------------------------------------------------------
*\

*Ramsey-reset test
estat ovtest
*Breeusch-Pagen Cook-Weisberg test for heteroskedasticity
estat hettest
*結果の保存
outreg2 using myreg.doc, replace ctitle(Model 1) 
est sto reg1

* Breusch Godfrey test 
estat bgodfrey, lags(1) 
matrix p1 = r(p) 
scalar rej1 = (p1[1,1]<alpha) 
estat dwatson 
scalar ddw = r(dw)

*----AICの比較から選択されない特定化
reg d_socgd  d_positive d_negative l.rnegative_socgd
*Akaike information criteria
estat ic

*\
Akaike's information criterion and Bayesian information criterion

-----------------------------------------------------------------------------
       Model |        Obs  ll(null)  ll(model)      df         AIC        BIC
-------------+---------------------------------------------------------------
           . |        656 -1511.357  -1482.033       4    2972.067   2990.012
-----------------------------------------------------------------------------
              
*\




*******フルモデル*******            

reg d_socgd  d_positive d_negative  d_cpi d_unemp d_nav d_iip  business_cycle crisis_dum l.rcpi_socgd
estat ic


*\
Akaike's information criterion and Bayesian information criterion

-----------------------------------------------------------------------------
       Model |        Obs  ll(null)  ll(model)      df         AIC        BIC
-------------+---------------------------------------------------------------
           . |        650 -1499.521  -1417.756      10    2855.512   2900.282
-----------------------------------------------------------------------------
*\


reg d_socgd  d_positive d_negative  d_cpi d_unemp d_nav d_iip  business_cycle crisis_dum l.runemp_socgd
estat ic


*\
Akaike's information criterion and Bayesian information criterion

-----------------------------------------------------------------------------
       Model |        Obs  ll(null)  ll(model)      df         AIC        BIC
-------------+---------------------------------------------------------------
           . |        650 -1499.521  -1415.833      10    2851.666   2896.436
-----------------------------------------------------------------------------
*\


reg d_socgd  d_positive d_negative  d_cpi d_unemp d_nav d_iip  business_cycle crisis_dum l.rnav_socgd
estat ic


*\
Akaike's information criterion and Bayesian information criterion

-----------------------------------------------------------------------------
       Model |        Obs  ll(null)  ll(model)      df         AIC        BIC
-------------+---------------------------------------------------------------
           . |        650 -1499.521  -1418.276      10    2856.552   2901.322
-----------------------------------------------------------------------------
*\


reg d_socgd  d_positive d_negative  d_cpi d_unemp d_nav d_iip  business_cycle crisis_dum l.riip_socgd
estat ic


*\
Akaike's information criterion and Bayesian information criterion

-----------------------------------------------------------------------------
       Model |        Obs  ll(null)  ll(model)      df         AIC        BIC
-------------+---------------------------------------------------------------
           . |        650 -1499.521  -1417.508      10    2855.017   2899.787
-----------------------------------------------------------------------------
*\


reg d_socgd  d_positive d_negative  d_cpi d_unemp d_nav d_iip  business_cycle crisis_dum l.rldp_socgd
estat ic


*\
Akaike's information criterion and Bayesian information criterion

-----------------------------------------------------------------------------
       Model |        Obs  ll(null)  ll(model)      df         AIC        BIC
-------------+---------------------------------------------------------------
           . |        650 -1499.521  -1415.999      10    2851.998   2896.767
-----------------------------------------------------------------------------
*\


*----AICの比較から選択される特定化
reg d_socgd  d_positive d_negative  d_cpi d_unemp d_nav d_iip  business_cycle crisis_dum l.rpositive_socgd
estat ic

*\
Akaike's information criterion and Bayesian information criterion

-----------------------------------------------------------------------------
       Model |        Obs  ll(null)  ll(model)      df         AIC        BIC
-------------+---------------------------------------------------------------
           . |        649 -1497.569  -1415.066      10    2850.133   2894.887
-----------------------------------------------------------------------------
*\

*Ramsey-reset test
estat ovtest*Ramsey-reset test
estat ovtest
*Breeusch-Pagen Cook-Weisberg test for heteroskedasticity
estat hettest
*結果の保存
outreg2 using myreg.doc, append ctitle(Model 2)
outreg2 using replication_sa.doc, replace ctitle(Poisitive)
est sto reg2

reg d_socgd  d_positive d_negative  d_cpi d_unemp d_nav d_iip  business_cycle crisis_dum l.rnegative_socgd
estat ic

*\
Akaike's information criterion and Bayesian information criterion

-----------------------------------------------------------------------------
       Model |        Obs  ll(null)  ll(model)      df         AIC        BIC
-------------+---------------------------------------------------------------
           . |        649 -1497.569  -1416.665      10    2853.331   2898.085
-----------------------------------------------------------------------------
*\

#------------------
# データの読み込み
#------------------
dat3_3<-read.csv("fecm_economy_result.csv", fileEncoding = "SJIS")


#----------------
# プロットの作成
#----------------
pd <-position_dodge(0.4)
df_full$モデル = factor(df_full$モデル, levels=c('フルモデル', '党派性無'))

pd <-position_dodge(0.4)
df_full$従属変数 = factor(df_full$従属変数, levels=c('社会志向(+)','社会志向(-)', '個人志向(+)',
                                             '個人志向(-)'))


pd <-position_dodge(0.4)
df_full$独立変数 = factor(df_full$独立変数, levels=c('否定情報', '肯定情報', '無党派割合', '野党支持率', '与党支持率'))

windows(40, 35)
ggplot(df_full, aes(x=独立変数, y=as.numeric(coef), color = モデル), na.rm=T) + 
  geom_errorbar(aes(ymin=as.numeric(lower), ymax=as.numeric(upper)), size=1, position=pd, width=.1, na.rm=T) +
  geom_line(position=pd, na.rm=T) +
  geom_point(position=pd, size=3, shape=21, fill="white") + 
  scale_colour_grey(name="モデル") +
  facet_wrap(~ 従属変数, ncol=2) +
  coord_flip() +
  geom_hline(yintercept = 0) +
  ylab("係数") + theme_bw() +
  geom_text(data = df_full[df_full$モデル == 'フルモデル', ], 
            aes(label = paste0(round(as.numeric(coef), 2), "\n[", round(as.numeric(lower), 2), ", ", round(as.numeric(upper), 2), "]")), 
            position = pd, vjust = 1.5, hjust = 0.5, size = 3) + 
  ggtitle("")

#----------------
# 関数の読み込み
#----------------
source("plot4_1.R")
source("plot4_2.R")
source("plot4_1anno.R")
source("plot4_2anno.R")


#------------------
# データの読み込み
#------------------
dat4_1<-read.csv("approve_economy.csv",fileEncoding = "SJIS")
attach(dat4_1)



#----------------------------
# 政党支持率のプロットの作成
#----------------------------
d1 <- plot4_1(ruling, "与党支持率", "与党支持率の推移と政治経済的出来事")
d2 <- plot4_1(dat4_1$opposite, "野党支持率", "野党支持率の推移と政治経済的出来事")
d3 <- plot4_1(indep, "無党派層割合", "無党派層割合の推移と政治経済的出来事")

windows(80, 50)
grid.arrange(d1, d2, d3, ncol = 1)


#--------------------------------------
# 内閣支持率・不支持率のプロットの作成
#--------------------------------------
d4 <- plot4_2(approve, "内閣支持率", "内閣支持率の推移と政治経済的出来事")
d5 <- plot4_2(disapprove, "内閣不支持率", "内閣不支持率の推移と政治経済的出来事")

windows(80, 60)
grid.arrange(d4, d5, ncol = 1)

#---------------------- 
# FIの計算：与党支持率
#---------------------- 
fdGPH(na.omit(ruling))

#---------------------- 
# FIの計算：野党支持率
#---------------------- 
fdGPH(na.omit(opposite))

#---------------------- 
# FIの計算：無党派割合
#---------------------- 
fdGPH(na.omit(indep))

#---------------------- 
# FIの計算：内閣支持率
#---------------------- 
fdGPH(na.omit(approve))

#------------------------ 
# FIの計算：内閣不支持率
#------------------------ 
fdGPH(na.omit(disapprove))

#------------------------
# 関数の読み込み
#------------------------
source("plot4_4df.R")
source("plot4_4_1.R")
source("plot4_4_2.R")

#------------------------
# データの読み込みと整形
#------------------------
dat4_4 <- read_csv("approve_economy.csv", locale = locale(encoding = "SJIS"))


#------------------- 
# 政権交代期の設定
#-------------------
admin <- data.frame(
  start = as.Date(c('1993-08-01', '2009-09-01')), 
  end   = as.Date(c('1993-05-01', '2012-11-01'))
)

#----------------------------
# データ作成とプロットの作成
#----------------------------
# 党派別の内閣支持率・不支持率の2系列の描画
tb_long_ap <- dat4_4 %>% with(plot4_4df(yearmon, ap_ldp, ap_opposite, "党派別の内閣支持率"))
d1 <- plot4_4_1(tb_long_ap, "与党派・野党派間の内閣支持率", "割合", "党派別の内閣支持率")

tb_long_dis <- dat4_4 %>% with(plot4_4df(yearmon, dis_ldp, dis_opposite, "党派別の内閣不支持率"))
d2 <- plot4_4_1(tb_long_dis, "与党派・野党派間の内閣不支持率", "割合", "党派別の内閣不支持率")

# 党派性差異の1系列の描画（支持率）
df_ap_diff <- dat4_4 %>% transmute(yearmon, 党派性差異 = abs(ap_ldp - ap_opposite))
d3 <- plot4_4_2(df_ap_diff, "内閣支持率の党派性差異", "党派性差異")

# 党派性差異の1系列の描画（不支持率）
df_dis_diff <- dat4_4 %>% transmute(yearmon, 党派性差異 = abs(dis_ldp - dis_opposite))
d4 <- plot4_4_2(df_dis_diff, "内閣不支持率の党派性差異", "党派性差異")

#---------------
# プロットの表示
#---------------
windows(70, 50)
grid.arrange(d1, d2, d3, d4, ncol = 1)

**************************
*  結果変数：与党支持率  *
**************************

***---1st step: 残差の予測値の計算
reg ruling sociotropic_gd
predict rsociogd_rul, resid

reg ruling liv_gd
predict rlivgd_rul, resid

reg ruling sociotropic_bd
predict rsociobd_rul, resid

reg ruling liv_bd
predict rlivbd_rul, resid

reg ruling cpi
predict rcpi_rul, resid

reg ruling unemp
predict runemp_rul, resid

reg ruling ln_nav
predict rnav_rul, resid

reg ruling ruling
predict rldp_rul, resid

reg ruling coincidentdi
predict rdi_rul, resid


reg ruling iip
predict riip_rul, resid

reg ruling interest
predict rint_rul, resid

reg ruling positive
predict rpositive_rul, resid

reg ruling negative
predict rnegative_rul, resid

reg ruling sociotropic_gd liv_gd positive negative  cpi unemp lnNikkei lnIIP  business_cycle crisis_dum corruption abuse admin_event  
predict rrul_gd, resid


reg ruling sociotropic_bd liv_bd positive negative  cpi unemp lnNikkei lnIIP  business_cycle crisis_dum corruption abuse admin_event  
predict rrul_bd, resid


***---2nd step: d値によるフィルタリングの実施
arfima d.rrul_gd, iterate(2000)
predict d_rrul_gd, fdifference


arfima d.rrul_bd, iterate(2000)
predict d_rrul_bd, fdifference



***---3rd step: OLSによるFECM推定
********
*図5－1*
********
*肯定的評価に関する採用されたモデル（選別の手順については図3-3コードを参照）
reg d_rul d_socgd d_livgd d_positive d_negative  d_cpi d_unemp d_nikkei d_iip  business_cycle crisis_dum corruption abuse admin_event  l.d_rrul_gd

*Akaike information criteria
estat ic
*Ramsey-reset test
estat ovtest
*Breeusch-Pagen Cook-Weisberg test for heteroskedasticity
estat hettest
*結果の保存
outreg2 using myreg.doc, replace ctitle(Model 1) 
est sto reg1

* Breusch Godfrey test 
estat bgodfrey, lags(1) 
matrix p1 = r(p) 
scalar rej1 = (p1[1,1]<alpha) 
estat dwatson 
scalar ddw = r(dw)

*否定的評価に関する採用されたモデル（同様に手順は図3－3を参照）
reg d_rul d_socbd d_livbd d_positive d_negative  d_cpi d_unemp d_nikkei d_iip  business_cycle crisis_dum corruption abuse admin_event  l.d_rrul_bd
***検定など省略


********
*図5－2*
********
*肯定的評価  
*1960-1980
reg d_rul d_socgd d_livgd d_positive d_negative  d_cpi d_unemp d_nikkei d_iip  business_cycle crisis_dum corruption abuse admin_event  l.d_rrul_gd if time_id >0 & time_id < 194
***検定など省略


*1980-2000
reg d_rul d_socgd d_livgd d_positive d_negative  d_cpi d_unemp d_nikkei d_iip  business_cycle crisis_dum corruption abuse admin_event  l.d_rrul_gd if time_id >193 & time_id < 435
***検定など省略

*2000-
reg d_rul d_socgd d_livgd d_positive d_negative  d_cpi d_unemp d_nikkei d_iip  business_cycle crisis_dum corruption abuse admin_event  l.d_rrul_gd if time_id >434
***検定など省略


*否定的評価
*1960-1980
reg d_rul d_socbd d_livbd d_positive d_negative  d_cpi d_unemp d_nikkei d_iip  business_cycle crisis_dum corruption abuse admin_event  l.d_rrul_bd if time_id >0 & time_id < 194
***検定など省略

*1980-2000
reg d_rul d_socbd d_livbd d_positive d_negative  d_cpi d_unemp d_nikkei d_iip  business_cycle crisis_dum corruption abuse admin_event  l.d_rrul_bd if time_id >193 & time_id < 435
***検定など省略

*2000-
reg d_rul d_socbd d_livbd d_positive d_negative  d_cpi d_unemp d_nikkei d_iip  business_cycle crisis_dum corruption abuse admin_event  l.d_rrul_bd if time_id >434
***検定など省略

****************
*  内閣支持率　*
****************

***---1st step: 残差の予測値の計算
reg approve sociotropic_gd
predict rsociogd_ap, resid

reg approve liv_gd
predict rlivgd_ap, resid

reg approve cpi
predict rcpi_ap, resid

reg approve unemp
predict runemp_ap, resid

gen ln_nav=log(nav_high)
reg approve ln_nav
predict rnav_ap, resid

reg approve ruling
predict rldp_ap, resid

reg approve coincidentdi
predict rdi_ap, resid

reg approve m2
predict rm2_ap, resid

reg approve iip
predict riip_ap, resid

reg approve interest
predict rint_ap, resid

reg approve positive
predict rpositive, resid

reg approve negative
predict rnegative, resid

reg approve ruling sociotropic_gd liv_gd positive negative  cpi unemp lnNikkei lnIIP  business_cycle crisis_dum corruption abuse admin_event  
predict rfull, resid


reg approve sociotropic_gd liv_gd positive negative  cpi unemp lnNikkei lnIIP  business_cycle crisis_dum corruption abuse admin_event  
predict rsocfull, resid


***---2nd step: d値によるフィルタリングの実施

arfima d.approve, iterate(2000)
predict d_ap, fdifference

arfima d.sociotropic_gd,  iterate(2000)
predict d_socgd, fdifference

arfima d.liv_gd, iterate(2000)
predict d_livgd, fdifference

arfima d.cpi, iterate(2000)
predict d_cpi, fdifference

arfima d.unemp, iterate(2000)
predict d_unemp, fdifference

arfima d.lnNikkei, iterate(2000)
predict d_nikkei, fdifference

arfima d.ruling, iterate(2000)
predict d_rul, fdifference

arfima d.lnIIP, iterate(2000)
predict d_iip, fdifference


arfima d.interest, iterate(2000)
predict d_int, fdifference

arfima positive, iterate(2000)
predict d_positive, fdifference

arfima negative, iterate(2000)
predict d_negative, fdifference

arfima d.indep, iterate(2000)
predict d_indep, fdifference

arfima rfull, iterate(2000)
predict d_rfull, fdifference

arfima rsocfull, iterate(2000)
predict d_rsocfull, fdifference

arfima rcsifull, iterate(2000)
predict d_rcsifull, fdifference


***---3rd step: OLSによるFECM推定
********
*図5－3*
********  
reg d_ap d_rul d_socgd d_livgd d_positive d_negative  d_cpi d_unemp d_nikkei d_iip  business_cycle crisis_dum corruption abuse admin_event  l.d_rfull
***検定など省略

********
*図5－4*
********
*1960-1980
reg d_ap d_rul d_socgd d_livgd d_positive d_negative  d_cpi d_unemp d_nav  business_cycle crisis_dum corruption abuse admin_event  l.d_rfull if time_id > 0 & time_id < 194
***検定など省略

*1980-2000
reg d_ap d_rul d_socgd d_livgd d_positive d_negative  d_cpi d_unemp d_nav  business_cycle crisis_dum corruption abuse admin_event  l.d_rfull if time_id > 193 & time_id < 435
***検定など省略

*2000-
reg d_ap d_rul d_socgd d_livgd d_positive d_negative  d_cpi d_unemp d_nav  business_cycle crisis_dum corruption abuse admin_event  l.d_rfull if time_id > 434 
est sto reg40
***検定など省略


*与党支持率なしモデルでは、上記から与党支持率を除いた特定化のもとに推定を行っている。内閣不支持率を結果変数とする各種モデルの手順は省略しているが、同様の手順により推定を行っている。