#!/usr/bin/env python # coding: utf-8 # # 移动平均(MA)模型 # ## 原理讲解 # ### 模型介绍 # # 另一类时间序列模型为“移动平均过程”（Moving Average Process，简记 MA）。记一阶移动平均过程为 MA（1）： # # $$y_{t}=\mu+\varepsilon_{t}+\theta \varepsilon_{t-1}\tag{1}$$ # # 其中，{$\varepsilon_{t}$} 为白噪声，而 $\varepsilon_{t}$ 的系数被标准化为1。由于 $y_t$ 可以被看成是白噪声的移动平均，故得名。考虑使用条件 MLE 估计。为此，假设 {$\varepsilon_{t}$} 为独立同分布且服从正态分布 $N\left(0, \sigma_{\varepsilon}^{2}\right)$。如果已经知道 {$\varepsilon_{t-1}$} ，则 # # $$y_{t} | \varepsilon_{t-1} \sim N\left(\mu+\theta \varepsilon_{t-1}, \sigma_{\varepsilon}^{2}\right)\tag{2}$$ # # 假设 $\varepsilon_0=0$，则可以知道 $\varepsilon_{1}=y_{1}-\mu$。给定 $\varepsilon_{1}$，则可以知道 $\varepsilon_{2}=y_{2}-\mu-\theta \varepsilon_{1}$。以此类推，则可以使用递推公式“$\varepsilon_{t}=y_{t}-\mu-\theta \varepsilon_{t-1}$”计算出全部 $\left\{\varepsilon_{1}, \varepsilon_{2}, \cdots, \varepsilon_{T}\right\}$。这样，在给定 $\varepsilon_0=0$ 的条件下，就可以写下样本数据 $\left\{y_{1}, y_{2}, \cdots, y_{T}\right\}$ 的似然函数，然后使用数值方法求解其最大化问题。 # # 更一般地，对于 $q$ 阶移动平均过程，记为 MA$(q)$： # # $y_{t}=\mu+\varepsilon_{t}+\theta_{1} \varepsilon_{t-1}+\theta_{2} \varepsilon_{t-2}+\cdots+\theta_{q} \varepsilon_{t-q}\tag{3}$ # # 也可以进行条件 MLE 估计，即在给定“$\varepsilon_{0}=\varepsilon_{-1}=\cdots=\varepsilon_{-q+1}=0$”的条件下，最大化样本数据的似然函数。 # ### MA(1) 模型最大似然函数推导 # # $y_t=μ+ε_t+θε_{t-1}$ # # $y_1=μ+ε_1+θε_0=μ+ε_1  ⇒ε_1=y_1-μ  ⇒y_1∼N(μ,σ_ε^2)$ # # $y_2=μ+ε_2+θε_1=(1-θ)μ+θy_1+ε_2  ⇒ε_2= y_2-θy_1-(1-θ)μ  ⇒y_2∼N((1-θ)μ+θy_1,σ_ε^2)$ # # 以此类推：$y_3∼N((1-θ+θ^2)μ+θy_2-θ^2 y_1,σ_ε^2)$ # # 因此： # # $$f\left(y_{1}\right)=\frac{1}{\sqrt{2 \pi \sigma_{\varepsilon}^{2}}} e^{-\frac{\left(y_{1}-\mu\right)^{2}}{2 \sigma_{\varepsilon}^{2}}}, f_{y_{2} | y_{1}}\left(y_{2} | y_{1}\right)=\frac{1}{\sqrt{2 \pi \sigma_{\varepsilon}^{2}}} e^{-\frac{\left(y_{2}-(1-\theta) \mu-\theta_{1}\right)^{2}}{2 \sigma_{\varepsilon}^{2}}}\tag{4}$$ # # $$\ln L=-\frac{1}{2} \ln 2 \pi-\frac{1}{2} \ln \sigma_{\varepsilon}^{2}-\frac{\left(y_{1}-\mu\right)^{2}}{2 \sigma_{\varepsilon}^{2}}-\frac{T-1}{2} \ln 2 \pi-\frac{T-1}{2} \ln \sigma_{\varepsilon}^{2}-\sum_{t=2}^{T} \frac{\left(y_{t}-c_{1} \mu-c_{2}\right)}{2 \sigma_{\varepsilon}^{2}}\tag{5}$$ # # 其中：$c_1=1-θ+θ^2-θ^3+......;    c_2=θy_(t-1)-θ^2 y_(t-2)+θ^3 y_(t-3)+......$ # ## statsmodels 库实现 # 详情请参考：[ARMA 模型](https://www.statsmodels.org/stable/generated/statsmodels.tsa.arima_model.ARMA.html#statsmodels.tsa.arima_model.ARMA) # In[1]: from statsmodels.tsa.arima_model import ARMA import pandas as pd data = pd.read_excel('../数据/上证指数与沪深300.xlsx') res = ARMA(data['sz'], order=(0,1)).fit() res.summary() # ## matlab实现 # # MAlnL.m # ```matlab # function lnL = MAlnL(Beta,Y) # %Beta(1) = u ;Beta(2) = sita;Beta(3) = sigma2; # T = length(Y); # lnL1 = -0.5*log(2*pi)-0.5*log(Beta(3)^2)-(Y(1)-Beta(1))^2/(2*Beta(3)^2)-0.5*log(2*pi)*(T-1)-0.5*(T-1)*log(Beta(3)^2); # for i = 2:T # for j = 1:i # c11(j) = (-Beta(2))^(j-1); # end # c1 = sum(c11)*Beta(1); # for j = 1:i-1 # c22(j) = Y(i-j)*(Beta(2)^j)*(-1)^(j+1); # end # c2 = sum(c22); # lnL2(i) = (Y(i)- c1 - c2)^2; # end # lnL = lnL1-0.5*sum(lnL2)/Beta(3)^2; # lnL = -lnL; # # % %导入数据 # % data = xlsread('../数据/上证指数与沪深300.xlsx'); # % f1 = data(:,1); f2 = data(:,2); # % # % % 2.初始参数设定 # % maxsize = 2000; % 生成均匀随机数的个数(用于赋初值) # % REP = 100; % 若发散则继续进行迭代的次数 # % nInitialVectors = [maxsize, 3]; % 生成随机数向量 # % MaxFunEvals = 5000; % 函数评价的最大次数 # % MaxIter = 5000; % 允许迭代的最大次数 # % options = optimset('LargeScale', 'off','HessUpdate', 'dfp','MaxFunEvals', ... # % MaxFunEvals, 'display', 'on', 'MaxIter', MaxIter, 'TolFun', 1e-6, 'TolX', 1e-6,'TolCon',10^-12); # % # % % 3.寻找最优初值 # % initialTargetVectors = unifrnd(0,10, nInitialVectors); # % RQfval = zeros(nInitialVectors(1), 1); # % for i = 1:nInitialVectors(1) # % RQfval(i) = MAlnL (initialTargetVectors(i,:), f1); # % end # % Results = [RQfval, initialTargetVectors]; # % SortedResults = sortrows(Results,1); # % BestInitialCond = SortedResults(1,2: size(Results,2)); # % # % % 4.迭代求出最优估计值Beta # % [Beta, fval exitflag] = fminsearch(' MAlnL ', BestInitialCond,options,f1); # % for it = 1:REP # % if exitflag == 1, break, end # % [Beta, fval exitflag] = fminsearch(' MAlnL ', BestInitialCond,options,f1); # % end # % if exitflag~=1, warning('警告：迭代并没有完成'), end # # ``` # ## Eviews 估计 # 估计方程 `y c ma(p)` #