文本分析与叙事冲击
Text Analysis & Narrative Shocks

# ==============================================================
# 文本分析与叙事冲击：完整可运行演示
# 依赖: pip install gensim scikit-learn jieba pandas numpy statsmodels
# ==============================================================
import re, numpy as np, pandas as pd
from gensim import corpora
from gensim.models import LdaModel
from sklearn.feature_extraction.text import TfidfVectorizer
import statsmodels.api as sm

# ---------- 0. 模拟央行会议纪要语料（真实研究替换为抓取文本） ----------
# 每条 = 一次会议纪要。我们故意混入"宽松派"与"紧缩派"措辞
np.random.seed(7)
dovish_words   = ["ease","accommodative","support","liquidity","cut","growth","downside","stimulus"]
hawkish_words  = ["tighten","inflation","overheating","rates","hawkish","anchor","overvalued","risk"]
neutral_words  = ["monitor","stable","expect","policy","committee","meeting"," outlook","data"]

def make_doc(tilt):
    # tilt=+1 偏鹰, -1 偏鸽, 0 中性
    bag = []
    for _ in range(60):
        u = np.random.rand()
        if tilt > 0 and u < 0.6:
            bag.append(np.random.choice(hawkish_words))
        elif tilt < 0 and u < 0.6:
            bag.append(np.random.choice(dovish_words))
        else:
            bag.append(np.random.choice(neutral_words))
    return " ".join(bag)

docs = [make_doc(np.random.choice([-1,0,1])) for _ in range(400)]

# ---------- 1. 预处理：分词 + 去停用词 ----------
stopwords = set(["the","a","an","of","to","and","is","in","for","on","our","we"])
def tokenize(text):
    toks = re.findall(r"[a-z]+", text.lower())
    return [t for t in toks if t not in stopwords and len(t) > 2]

tokenized = [tokenize(d) for d in docs]
dictionary = corpora.Dictionary(tokenized)
dictionary.filter_extremes(no_below=5, no_above=0.5)   # 过滤极端词
corpus_bow = [dictionary.doc2bow(t) for t in tokenized]

# ---------- 2. TF-IDF 加权 ----------
tfidf = TfidfVectorizer(tokenizer=lambda x: x, lowercase=False,
                        token_pattern=None)
tfidf_matrix = tfidf.fit_transform([" ".join(t) for t in tokenized])
print("文档-词矩阵维度:", tfidf_matrix.shape)

# ---------- 3. LDA 主题模型（gensim） ----------
K = 3   # 主题数；实际用 coherence 选择
lda = LdaModel(corpus=corpus_bow, id2word=dictionary,
               num_topics=K, random_state=42,
               alpha='auto', eta='auto',
               passes=20, iterations=200)

print("\n=== 各主题 top-8 词 ===")
for k in range(K):
    top = lda.show_topic(k, topn=8)
    print(f"主题{k}: " + ", ".join([w for w,_ in top]))

# 提取每篇文档的文档-主题分布 theta_d（后验推断）
theta = np.array([[p for _,p in lda.get_document_topics(bow, minimum_probability=0.0)]
                  for bow in corpus_bow])
print("\n文档-主题分布 theta 维度:", theta.shape)   # (D, K)

# ---------- 4. 词典法情绪指数（鸽 vs 鹰） ----------
hawk_set   = set(hawkish_words)
dove_set   = set(dovish_words)
def sentiment(tokens):
    n_pos = sum(t in hawk_set for t in tokens)
    n_neg = sum(t in dove_set for t in tokens)
    return (n_pos - n_neg) / max(len(tokens),1)

tone = np.array([sentiment(t) for t in tokenized])   # >0 偏鹰, <0 偏鸽
print("\n情绪指数 tone 描述: mean=%.3f std=%.3f" % (tone.mean(), tone.std()))

# ---------- 5. 把文本指数与利率做回归 ----------
# 模拟政策利率（真实研究用 shibor/DR007/Fed Funds）
policy_rate = 0.02 + 0.30*tone + 0.1*np.random.randn(len(tone))*0.001
df = pd.DataFrame({"tone":tone, "rate":policy_rate,
                   "theta_k0":theta[:,0], "theta_k1":theta[:,1]})
X = sm.add_constant(df[["tone","theta_k0","theta_k1"]])
ols = sm.OLS(df["rate"], X).fit()
print("\n=== 文本语气对政策利率的回归 ===")
print(ols.summary().tables[1])
# 解读: tone 系数应显著为正——偏鹰的纪要对应更高政策利率

# ==============================================================
# R: 用 stm 包估计央行纪要 LDA / 结构主题模型
# install.packages(c("stm","topicmodels","quanteda"))
# ==============================================================
library(quanteda); library(stm); library(topicmodels)

# ---------- 0. 读入语料：一行一篇会议纪要 ----------
# df <- read.csv("pboc_minutes.csv", stringsAsFactors=FALSE)
# df$date <- as.Date(df$date)
# 这里用内置数据集演示
data(gadarian)            # quanteda 内置：妇女参政论者传单文本
docs  <- gadarian$documents
meta  <- gadarian[ , c("treatment","pid","year","age","female","ideology")]

# ---------- 1. 构造文档-词矩阵（quanteda） ----------
corp  <- corpus(docs)
toks  <- tokens(corp, remove_punct=TRUE, remove_numbers=TRUE)
toks  <- tokens_tolower(toks)
toks  <- tokens_remove(toks, stopwords("en"))
toks  <- tokens_wordstem(toks, language="english")
dfm_  <- dfm(toks)
dfm_  <- dfm_trim(dfm_, min_termfreq = 5, docfreq_max = 0.5)

# ---------- 2. stm：把协变量放进主题先验 ----------
# prevalence = ~ treatment + ideology + s(age)
# 表示"主题占比"受 treatment、意识形态、年龄影响
stm_fit <- stm(documents  = dfm_,
               K          = 5,
               prevalence = ~ treatment + ideology,
               data       = meta,
               seed       = 2024)

# 看每个主题的代表词
labelTopics(stm_fit, n = 8)

# ---------- 3. 主题占比随协变量变化（宏观研究常用） ----------
# estimateEffect: 在控制其他变量后, treatment 对各主题占比的影响
eff <- estimateEffect(1:5 ~ treatment + ideology,
                      stm_fit, meta, uncertainty = "Global")
plot(eff, "treatment", cov.value1 = "treatment", cov.value2 = "control")

# ---------- 4. 输出文档-主题分布 theta 供后续回归 ----------
theta_hat <- stm_fit$theta          # D x K 矩阵
write.csv(theta_hat, "doc_topic_proportions.csv")

# ---------- 5. 备选：经典 LDA（topicmodels 包） ----------
# dtm <- convert(dfm_, to = "topicmodels")
# lda_fit <- LDA(dtm, k = 5, method = "Gibbs",
#                control = list(seed=123, burnin=1000, thin=100, iter=5000))
# posterior(lda_fit)$terms     # phi_k
# posterior(lda_fit)$topics    # theta_d