此文內容為《動態(tài)隨機一般均衡(DSGE)模型》的筆記,李向陽老師著,清華大學出版社出版。
我只將個人會用到的知識作了筆記,并對教材較難理解的部分做了進一步闡述。為了更易于理解,我還對教材上的一些部分(包括代碼和正文)做了修改。
僅供學習參考,請勿轉載,侵刪!
上一篇有關Dyanre的文章為大家介紹了有關確定性模擬的技術。在本次推文,我將介紹隨機模擬技術。另外,我還會介紹有關參數估計的問題。
3.9 隨機模擬分析:stock_simul
隨機模擬是Dyanre最常用的計算命令之一。隨機模擬首先使用擾動算法完成模型的求解任務,然后再求解的基礎上計算內生變量的脈沖響應、各階矩,最后進行隨機模擬,——即根據隨機抽取的外生沖擊的樣本,模擬出內生變量的增長路徑。
3.9.1 隨機模擬命令簡介
在Dynare中,隨機模擬的選項多達44個,常用的有:
drop = xx,計算內生變量的各階矩時,需要首先指定丟棄的期數,默認為100hp_filter = xx,計算內生變量的各階矩時,首先使用HP濾波,xx代表值
irf=xx,計算xx期的脈沖響應,默認為40irf_shocks=(shock1, shock2,...),針對給定的外生沖擊計算脈沖響應,默認計算全部irf_plot_threshold=xx,指定畫脈沖響應圖形的門檻大小,默認nocorr,指定不要在Matlab窗口顯示相關系數矩陣nofunctions,指定不要在Matlab窗口顯示政策函數和轉換方程nomoments,指定不要在Matlab命令窗口顯示各階矩nograph,指定不要顯示和保存圖形nodisplay,指定不要顯示圖形,但保存到硬盤noprint,指定不要顯示任何計算信息,對循環(huán)很有幫助order=xx,指定Taylor近似對階數,可為1``2``3,默認為2periods=xx,指定隨機模擬的期數,同時也是決定計算理論矩還是模擬矩的關鍵所在pruning,指定剪枝算法,即在迭代過程中丟棄高階項
3.9.2 隨機模擬示例
以第三章第一節(jié)的RBC模型為例,這里示例如何使用模擬結果中的水平變量進行做圖分析:
/*
* 《動態(tài)隨機一般均衡(DSGE)模型》教材
* 隨機模擬示例
* 2020年3月30日
* @愛吃漢堡薯條
* @侵刪
*/
%----- 源代碼9 -----%
var c k lab z I; % 內生變量聲明
varexo e; % 外生變量聲明
parameters beta theta delta alpha tau rho s; % 參數聲明
% 參數賦值
beta = 0.987;
theta = 0.357;
delta = 0.012;
alpha = 0.4;
tau = 2;
rho = 0.95;
s = 0.01;
%----- 源代碼10 -----%
model;
% 定義模型的“局部變量”,以提高編程的便利性
# mar_c = (c^theta * (1-lab)^(1-theta) ) ^ (1-tau); % 局部變量的聲明使用“#”命令
# mar_c1 = (c(+1)^theta * (1-lab(+1))^(1-theta) ) ^ (1-tau);
% (1) 歐拉方程
[name = 'Euler equation'] % 這是系統(tǒng)標簽(tag),可以對每個均衡條件進行標識
mar_c/c = beta * (mar_c1/c(+1)) * (1 + alpha*exp(z)*k^(alpha-1)*lab(+1)^(1-alpha) - delta);
% (2) 勞動供給方程
[name = 'Labor supply']
((1-theta)/theta) * (c/(1-lab)) = (1-alpha)*exp(z)*k(-1)^alpha*lab^(-alpha);
% (3) 資源約束方程
[name = 'Resource constraint']
c + k - (1-delta)*k(-1) = exp(z)*k(-1)^alpha*lab^(1-alpha);
% (4) 資本積累方程
[name = 'Capital accumulation']
I = k - (1-delta)*k(-1);
% (5) 外生沖擊AR(1)過程
[name = 'Technology shock']
z = rho*z(-1) + s*e;
end;
%----- 源代碼11 -----%
initval; % 初值模塊,設定內生變量和外生變量的初始值
k = 0.5;
c = 0.5;
lab = 0.5;
z = 0;
I = 0.5;
e = 0;
end;
shocks; % 外生沖擊設定,直接設定外生沖擊的方差
var e = s^2;
end;
steady; % 計算穩(wěn)態(tài)
resid; % 給定內生變量的取值(通常為穩(wěn)態(tài)),計算均衡條件對應的靜態(tài)方程的殘差
%----- 源代碼12 -----%
% 如果不指定模擬的次數“periods”,Dynare就不會進行模擬
stoch_simul(periods=1000, irf=400, order=1);
% 將模擬結果保存為“simudata.mat”文件
dynasave('simudata.mat');
則有下圖:

在這里,隨機模擬命令指定了3個選項:
模擬1000期
計算40期脈沖響應圖
進行二階泰勒近似求解
前面已經提到,內生變量的模擬路徑結果儲存于oo_.endo_simul結構數組中。其行的順序對應著以聲明順序排序的內生變量,其列則為模擬期數。
其中,技術沖擊的均值為:
>> mean(oo_.endo_simul(4,:))
ans = -7.2640e-05
3.10 參數估計簡介
模型中結構參數的賦值方法有兩種,分別是:
- 校準(Calibration),直接賦值
- 估計(Estimation),使用實際統(tǒng)計或觀測數據進行估計
簡單來說,Dynare對參數估計就是使用某些方法,比如貝葉斯估計、極大似然估計,從數據中發(fā)現真實取值
3.10.1 極大似然估計與貝葉斯估計的基本邏輯
本文僅簡單闡述極大似然估計和貝葉斯估計的基本原理和邏輯,對背后的復雜數理知識,請參考教材和博客。
3.10.1.1 極大似然估計
極大似然估計(Maximum Likelihood Estimation, MLE)就是選擇參數 來最大化如下的似然函數(Likelihood),其原理比較簡單:
其中,條件概率 是使用
濾波 進行計算的。
總而言之,它是選擇 使已經發(fā)生時間的概率最大。
3.10.1.2 貝葉斯估計
貝葉斯估計(Bayesian Estimation)的基本原理非常直觀,就是在參數先驗分布(Prior)的基礎上,結合數據的信息(),找到參數的后驗分布(Posterior)。也就是說,貝葉斯估計是一種簡單的
映射,依據數據信息,將先驗分布映射為后驗分布。
由于先驗分布是依據經驗而假設的分布,可能存在一定的偏誤,而數據的作用其實是將數據中有用的信息融入,以糾正偏誤并形成后驗分布。

從邏輯上說,貝葉斯估計的背后就是貝葉斯法則,即參數的后驗分布由條件概率公式:
推導出:
從數據使用和先驗假設的角度看,貝葉斯估計位于極大似然估計和參數校準之間:
- 極大似然估計:直接使用數據中的信息,沒有先驗
- 貝葉斯估計:即使用數據中的信息,也具有先驗信息
- 校準:僅使用先驗信息,沒有使用數據的信息
貝葉斯估計的先驗分布是一般意義上的分布,即標準差不為0,不是一個點,而是一個區(qū)域。嚴格地說,是整個參數空間的一個子區(qū)域,作為候選區(qū)域。先驗分布可以理解為施加于似然函數的“權重”或者“偏好”,使得參數的選擇范圍集中于均值附近。
換句話說,先驗分布使得“巨大的參數選擇空間”變小,使似然函數集中于“更有意義”的子區(qū)域去尋找合適的參數估計值。Dynare的貝葉斯分布可分為兩步:
- 使用數值求解算法,計算似然函數取最大值的
眾數(Model) - 從
眾數開始,使用MCMC方法模擬后驗分布,然后計算各種感興趣的矩(Moments)
3.10.2 馬爾可夫鏈——蒙特卡洛(MCMC)算法
由于參數的后驗分布(密度函數)在大多數DSGE模型中(甚至某些RBC模型中)無法使用解析式的形式來表達,因此不能通過貝葉斯法則計算出來,只能通過隨機抽樣的方法,即蒙特卡洛算法(Monte Carlo)來近似逼近后驗分布。
馬爾可夫鏈-蒙特卡洛方法(Markov Chain-Monte Carlo, MCMC)是一類方法的統(tǒng)稱,是貝葉斯估計的邏輯方法和框架,非具體估計的算法。可以參考 知乎 網友的文章。實際上我們不太需要了解背后的數學原理。
關于本節(jié),我們最重要的是基本了解Dynare貝葉斯估計的核心邏輯,如圖3.12所示:

3.10.3 Dynare參數估計和一個例子
我們著重于介紹如何用Dynare內置的估計命令,對一個或多個參數進行編程估計。
3.10.3.1 Dynare參數估計的基本邏輯
Dynare參數估計是使用以下模型的一階線性化系統(tǒng)進行估計:
其中, 是待估計參數向量;
是外生沖擊;
是內生變量組成的向量。
Dynare求解參數估計的基本邏輯是:
- 計算內生變量
的穩(wěn)態(tài)值
- 線性化系統(tǒng)均衡條件
- 求解線性化模型
- 使用Kalman濾波算法計算對數似然函數(Log-Likelihood)
- 找到極大似然值分布的眾數(Model)
- 使用Metropolis-Hastings算法模擬后驗分布
3.10.3.2 Dynare參數估計的語法
Dyanre參數估計的基本語法分成三部分:
- 聲明參數估計模塊
- 待估參數本身的聲明
- 使用
estimation命令
Dynare中聲明估計參數模塊的內置命令是以estimated_params;開頭,以end;結束,中間輸入待估計參數聲明。待估計參數有兩種常見的類型:
- 結構參數
- 外生參數的標準差參數
具體實例如下:
estimated_params; % 結構參數的估計
gamma, normal_pdf, 1,0.05;
alpha, uniform_pdf, , , 0, 1; % uniform分布需要留逗號,見下文
end;
estimated_params; % 外生沖擊標準差參數的估計
stderr e, inv_gamma_pdf, 0.01, inf;
end;
結構參數的聲明是:
- 以待估計參數的名稱開頭
- 然后是該參數服從的先驗分布的名稱,該名稱由Dynare預先定義,如表3.11所示:
| 先驗分布名稱 | 符號形式 | 分布取值范圍 |
|---|---|---|
beta_pdf |
||
gamma_pdf |
||
inv_gamma_pdf |
||
normal_pdf |
||
uniform_pdf |
其中, 和
分別表示均值和標準差;參數
、
僅用于
Beta分布、 Gamma分布和Uniform分布,可以為無窮大。此外,Uniform分布因為不需要聲明均值和標準差,因此需要在參數聲明中留出兩個逗號隔開的空格,然后聲明 、
參數。
標準差參數的聲明:
- 需要使用關鍵字
stderr后跟外生沖擊的名稱,而不是外生沖擊對應的標準差的名稱 - 然后是指定其服從的先驗分布,聲明語法與前面的相同
3.10.3.2 注意事項
對于先驗分布的選擇,確實是一個困難的問題。選擇先驗分布具有主觀性,這成為貝葉斯分布受質疑和批評的主要原因。一般參考:
- 對于標準差
的估計,通常選擇
Inverse Gamma 分布 - 對于
過程的滯后參數
通常選擇
Beta分布 - 其他參數,大多使用
Normal分布
在估計之前,還要注意以下幾個問題:
數據問題:估計參數的前提是模型中某些內生變量是可以觀測到或者由觀測值經過處理得到,且觀測變量的個數不能超過模型中外生沖擊的個數,否則無法計算模型的似然函數。
觀測數據與內生變量的匹配:在估計時,要求模型文件的變量形式和觀測數據相匹配。
-
觀測變量的聲明問題:內生變量如果聲明為可觀測變量,則需要使用以下命令
% 語法 varobs [VARIABLE_NAME]; % 例如 varobs pi tau a; % 聲明三個內生變量為觀測變量
3.10.3.3 estimation命令介紹
Dynare會在estimation命令中調取外部數據(使用datafile選項指定外部數據文件),并加載到內存中去。觀測變量聲明命令varobs會告訴Dynare去尋找變量名為pi、tau、a的三個序列,如果找不到,會提示錯誤。
Estimation`命令及其選項問題:這里僅介紹幾個最重要的關鍵選項
-
datafile=xxx:可選擇Matlab的.m文件、.mat矩陣或者.csv文件等 -
conf_sig=xxx:定義置信區(qū)間的寬度(如xxx=.95定義了“95%的置信區(qū)間”) -
frst_obs=xxx:定義由load_Gali_est_data.m加載的數據,第一個開始用于對估計的數據位置,xxx必須是一個正整數,缺省為1 -
nobs=yy:Dynare估計中使用的觀測變量序列長度(t=xxx到t=yy+xxx-1,同樣的yy必須是正整數,缺省為全部觀測變量 -
mode_check:畫出后驗分布眾數附近的后驗分布密度圖 -
mode_compute=xxx:選擇優(yōu)化算法,取值0~10,缺省使用4(該命令在不斷更新中) -
forecast=n:在觀測變量后一期之后預測期值
-
smmother:計算內生變量和外生沖擊的后驗分布的相關矩 -
mh_replic=xxx:定義計算后驗分布時MCMC重復的次數,缺省為2000 -
mh_jscale=xxx:定義MCMC中的乘法因子。理想情況下應該選擇使“接受率”落在25%~33%之間。該參數需要多次調試后獲得比較合適的值,缺省為0.2 -
mh_nblocks=xxx:定義Metropolis-Hastings算法平行進行的馬爾可夫鏈的個數,缺省為2 -
filtered_vars:計算濾波變量,結果存儲在oo_.FilteredVarialbes中 -
moments_varendo:計算內生變量的后驗分布的理論矩,結果存儲在oo_.PosteriorTheoreticalMoments中
3.10.3.4 極大似然估計的例子
考慮一個新凱恩斯模型,不含投資但含粘性價格和產出缺口。
(1)模型
其中:
-
為消費跨期替代彈性的倒數
-
為勞動負效用的外生沖擊,服從
過程:$\tau_{t}=\lambda \tau_{t-1}+\epsilon_{t}^{\tau}, \epsilon_{t}^{\tau} \sim i.i.Dá
-
為勞動供給的
彈性的倒數
首先是新凱恩斯IS曲線(NKIS)和新凱恩斯菲利普斯曲線(NKPC):
其次是自然水平產出和自然利率方程:
假設貨幣沖擊政策為Taylor規(guī)則:
從而,模型由 和六個均衡條件:
- NKIS曲線
- NKPC曲線
- 自然利率方程
- Taylor規(guī)則
- 技術沖擊的
過程
- 偏好沖擊的
過程
接下來,對上述模型文件進行隨機模擬,模擬5000個樣本,作為下一步估計的數據。然后使用上述模擬的樣本作為觀測數據,對模型中3個結構參數進行估計,如表3.12所示,以示例如何在Dynare中分別使用極大似然估計和貝葉斯估計估計此三個參數:
| 待估計參數 | 含義 |
|---|---|
| 技術沖擊 |
|
| 勞動負效用外生沖擊 |
|
| 技術沖擊 |
隨機模擬的代碼如下(請首先運行):
// 生成用于參數估計的隨你模擬序列
//Written by Xiangyang Li@SCC
//declaration of the endogenous variables, log-devation form
var a // technology
pi //CPI inflation
i //nominal rate
istar //nature real rate
tau //labor disutility shock
x //output gap = output - natural output
dy //growth rate of output
;
varexo eps_a eps_tau eps_i;
parameters
sigma phi beta kappa phi_x phi_pi rho lambda rhoi theta;
// Parameter Values
sigma = 1; //inverse of inter-temporal elasticity of substitution of consumption
beta = 0.97; //discount factor of households
phi_x = .15; //coefficient of output gap in monetary policy rule
phi_pi = 1.5; //coefficient of inflation in monetary policy rule
rhoi = 0.8; //persistence of nominal interest rate;
rho = 0.8; //persistence of technology
lambda = 0.5; //persistence of labor disutility shock;
phi = 1; //inverse of Frisch elasticity of labor supply
theta = 0.75; //Calvo stickiness parameter
kappa = ((1-theta)*(1-beta*theta))/theta; //simplifying parameter
// The model in log-linear
model(linear);
//(1) Philiphs Curve - Calvo Pricing Equation
beta*pi(+1) + kappa*x = pi;
//(2) New Keynesian IS curve
sigma*(i - pi(+1)-istar) = x(+1) - x;
//(3) natural real rate definition
istar = sigma*(1+phi)*(rho-1)/(sigma+phi)*a + sigma*(1-lambda)/(sigma+phi)*tau;
//(4)Taylor Rule
i= rhoi*i(-1) + (1-rhoi)*(phi_pi*pi + phi_x*x) +eps_i;
//(5)Technology shock
a = rho*a(-1) + eps_a;
//(6) Labor disutility shock
tau = lambda*tau(-1) + eps_tau;
//(7) Growth rate of output - by output gap
dy = x - x(-1) + (1+phi)/(sigma+phi)*(a-a(-1)) - (tau - tau(-1))/(sigma+phi); // construction of output growth
end;
shocks;
var eps_a;
stderr .02;
var eps_tau;
stderr .02;
var eps_i;
stderr 0.02;
end;
set_dynare_seed=1;
stoch_simul(periods=5000, irf=7, nograph);
gcsim = oo_.endo_simul;
save data gcsim;
(2)最大似然估計的實踐
在下面的代碼中,使用了estimated_param;設定了初始值,使用estimated_param_bounds;設定了上下邊界,此外使用varobs定義了3個觀測變量(請在運行3.10.3.4的(1)后運行):
// Maximum Likelihood estimation of the Gali-Christiano Model
//we are going to estimate two persistence parameters, rho and lambda
// and one standard devation parameters
//Written by Xiangyang Li@SCC@2016.12
//Modified by Xiangyang Li@2017.10.12
//declaration of the endogenous variables, log-devation form
var a // technology
pi //CPI inflation
i //nominal rate
istar //nature real rate
tau //labor disutility shock
x //output gap = output - natural output
dy //growth rate of output
;
varexo eps_a eps_tau eps_i;
parameters sigma phi beta kappa phi_x phi_pi rhoi theta;
parameters rho lambda; //to be estimated
// Parameter Values
sigma = 1; //inverse of inter-temporal elasticity of substitution of consumption
beta = 0.99; //discount factor of households
phi_x = .15; //coefficient of output gap in monetary policy rule
phi_pi = 1.5; //coefficient of inflation in monetary policy rule
rhoi = 0.8; //persistence of nominal interest rate;
//rho = 0.8; //persistence of technology, to be estimated
//lambda = 0.5; //persistence of labor disutility shock,to be estimated;
phi = 1; //inverse of Frisch elasticity of labor supply
theta = 0.75; //Calvo stickiness parameter
kappa = ((1-theta)*(1-beta*theta))/theta; //simplifying parameter
// The model in log-linear
model(linear);
//(1) Philiphs Curve - Calvo Pricing Equation
beta*pi(+1) + kappa*x = pi;
//(2) New Keynesian IS curve
sigma*(i - pi(+1)-istar) = x(+1) - x;
//(3) natural real rate definition
istar = sigma*(1+phi)*(rho-1)/(sigma+phi)*a + sigma*(1-lambda)/(sigma+phi)*tau;
//(4)Taylor Rule
i= rhoi*i(-1) + (1-rhoi)*(phi_pi*pi + phi_x*x) +eps_i;
//(5)Technology shock
a = rho*a(-1) + eps_a;
//(6) Labor disutility shock
tau = lambda*tau(-1) + eps_tau;
//(7) Growth rate of output - by output gap
dy = x - x(-1) + (1+phi)/(sigma+phi)*(a-a(-1)) - (tau - tau(-1))/(sigma+phi); // construction of output growth
end;
// 設置沖擊
shocks;
var eps_tau; stderr 0.01;
var eps_i; stderr 0.01;
end;
// 極大似然估計fa:
// 初始化參數值,用于優(yōu)化
estimated_params;
stderr eps_a, 0.01; // 外生沖擊參數的估計聲明
rho, .80; // 結構參數的估計聲明
lambda, .50;
end;
// 設置參數的邊界
estimated_params_bounds;
stderr eps_a, 0.001, .2;
rho, .001,.95;
lambda, .001,.95;
end;
// 用于估計參數的觀測源:觀測變量
varobs pi tau a;
// 開始估計
estimation(datafile=load_Gali_est_data,conf_sig =.95,first_obs=101,forecast =46,nobs=4000,mode_check,mode_compute=10) a pi i istar tau x dy;
運行上面的代碼,意味著:
- 使用95%置信區(qū)間
- 從第101個數據開始作為觀測值
- 在觀測變量以后預測46期
- 使用的序列長度是4000期
- 畫出后驗分布眾數附近的后驗分布密度圖
- 使用第10種優(yōu)化算法(見Manual)
(3)關于觀測變量的分析
在默認情況下,極大似然估計的結果令人滿意,非常接近真實值;但選擇不同的觀測變量,會影響參數的點估計,如表3.13所示:

直觀上,選擇觀測數據時,要盡量選擇和待估計參數關系緊密的數據序列,包括使用盡可能多的數據項,以獲得準確的信息。
(4)關于序列長度的分析
同樣,觀測變量序列的長度也會影響估計。可以看出,序列長度縮小100倍時,標準差約增大十倍。除了參數 以外,另外兩個參數的估計值變化很小,說明極大似然估計具有一定的穩(wěn)健型:

3.10.3.5 貝葉斯估計的例子
使用貝葉斯分布,要先指定先驗分布。這里設定:
- 技術標準差參
數的均值和真實值相同,標準差為10
-
和
兩個參數的均值為真實值,標準差都為0.02
- 還設定了3個參數的上下界和初始值
如下所示(請在運行3.10.3.4的(1)后運行):
// Maximum Likelihood estimation of the Gali-Christiano Model
//we are going to estimate two persistence parameters, rho and lambda
// and one standard devation parameters
//Written by Xiangyang Li@SCC@2016.12
//Modified by Xiangyang Li@2017.10.12
//declaration of the endogenous variables, log-devation form
var a // technology
pi //CPI inflation
i //nominal rate
istar //nature real rate
tau //labor disutility shock
x //output gap = output - natural output
dy //growth rate of output
;
varexo eps_a eps_tau eps_i;
parameters sigma phi beta kappa phi_x phi_pi rhoi theta;
parameters rho lambda; //to be estimated
// Parameter Values
sigma = 1; //inverse of inter-temporal elasticity of substitution of consumption
beta = 0.99; //discount factor of households
phi_x = .15; //coefficient of output gap in monetary policy rule
phi_pi = 1.5; //coefficient of inflation in monetary policy rule
rhoi = 0.8; //persistence of nominal interest rate;
//rho = 0.8; //persistence of technology, to be estimated
//lambda = 0.5; //persistence of labor disutility shock,to be estimated;
phi = 1; //inverse of Frisch elasticity of labor supply
theta = 0.75; //Calvo stickiness parameter
kappa = ((1-theta)*(1-beta*theta))/theta; //simplifying parameter
// The model in log-linear
model(linear);
//(1) Philiphs Curve - Calvo Pricing Equation
beta*pi(+1) + kappa*x = pi;
//(2) New Keynesian IS curve
sigma*(i - pi(+1)-istar) = x(+1) - x;
//(3) natural real rate definition
istar = sigma*(1+phi)*(rho-1)/(sigma+phi)*a + sigma*(1-lambda)/(sigma+phi)*tau;
//(4)Taylor Rule
i= rhoi*i(-1) + (1-rhoi)*(phi_pi*pi + phi_x*x) +eps_i;
//(5)Technology shock
a = rho*a(-1) + eps_a;
//(6) Labor disutility shock
tau = lambda*tau(-1) + eps_tau;
//(7) Growth rate of output - by output gap
dy = x - x(-1) + (1+phi)/(sigma+phi)*(a-a(-1)) - (tau - tau(-1))/(sigma+phi); // construction of output growth
end;
shocks;
var eps_tau;
stderr 0.01;
var eps_i;
stderr 0.01;
end;
// 貝葉斯估計
// 設置先驗分布
estimated_params;
// eps_a的均值為0.02,標準差為10
stderr eps_a, inv_gamma_pdf, 0.02, 10;
// 設置另外兩個參數的先驗
rho, beta_pdf, 0.80, 0.04;
lambda, beta_pdf, 0.80, 0.04;
end;
estimated_params_bounds;
stderr eps_a, 0.001, .2;
rho, .001,.95;
lambda, .001,.95;
end;
// 設置初值,MLE則不需要這一模塊
estimated_params_init;
stderr eps_a, 0.019;
rho, 0.82;
lambda, 0.51;
end;
// 設置觀測變量
varobs pi tau a;
estimation(datafile=load_Gali_est_data, // 數據文件
conf_sig =.95, // 95%置信區(qū)間
first_obs=101, // 從101位開始使用數據
forecast =46, // 在觀測序列后預測46期
nobs=40, // 使用40期數據進行預測
mode_check, // 畫出后驗分布分布密度圖
mode_compute=4, // 使用第4種優(yōu)化算法
mh_replic=1200, // MCMC重復次數
mh_jscale=1.4, // MCMC的乘法因子
mh_nblocks=2) // 馬爾可夫鏈個數
a pi i istar tau x dy;
運行上面的代碼,可以得到幾個圖,其中下圖:

可以直觀看出貝葉斯估計的先驗和后驗的情況。
貝葉斯估計比極大似然估計更加穩(wěn)健(請看教材),所以一般使用貝葉斯估計。
ng$