网站备案 拨测百度搜索引擎seo

导读：一般人都很喜欢用Stata来构建面板数据模型，我Stata相对较弱，更喜欢Matlab，这里就给出用Matlab来进行面板数据的混合效应建模实例，特别是code解读

面板数据的四大好处

面板数据(Panel Data)能够从时间和截面构成的二维空间来反映数据的变化规律，具有控 制个体的异质性、减少回归变量之间的多重共线性等优点，从而开始被广泛地应用于经济研究中，成为目前计量经济学领域研究的热点问题之一。

面板数据是指对总体中的给定样本在一段特定时期进行多重观察的数据集。这种多重观察，既包括对样本单位在某一时期 (时点)上多个特性的观察，也包括对该样本单位的这些特性在一段时期内的连续观察。面板数据早期主要用于天文和农业等领域，在1932年由Mundlak、Baleslfa和Nerlovc引入到计量经济领域。

参见 MarcNerlove.  An essay on the history ofPanel Data . Econometrics ,2000

与此同时，面板数据方法的研究亦逐渐成为目前计量经济模型研究的热点之一，许多专家和学者都提出了自己的看法及相应的检验方法，这主要是由于面板数据能够从时间和截面构成的二维空间来反映数据的变化规律，从而它与纯时间序列和截面数据相比有许多优点，归纳起来有以下几点：

(1) 面板数据能够控制个体的异质性。面板数据库显示个体(包括个人、企业、地区或国家)之间存在的差异，而单独的时间序列和横截面数据不能有效地反映这种差异。如果只是简单地使用时间序列和横截面数据进行分析，就可能得到有偏的结果。此外，面板数据分析熊够控制在时间序列和横截面研究中不能控制的涉及地区和时间为常数的情况，也就是说当个体在时间或地区分布中存在着常数的变量(例如受教育程度、电视广告等)时，如果在模型中不考虑这些变量，有可能会的有偏的结果。

(2) 面板数据能够提供更多信息，减少回归变量之间的多重共线性，增加自由度，从而提高参数估计的有效性。

(3)面板数据能够更好地研究动态调节。使用横截面数据进行分析，看上去相对稳定但却隐藏了许多变化，面板数据分析由于包含较长时间，能够进行前后对比，从而可以弄清如经济政策变化对我国经济状况的影响等问题。

(4)相对于纯时间序列与截面数据而言，使用面板数据能够构造和检验更复杂的行为模型；此外，面板数据可以收集到更准确的微观单位(个人、企业、家庭)的信息，由此得到的总体数据可以消去测量误差的影响。

一般人都很喜欢用Stata来构建面板数据模型，我Stata相对较弱，更喜欢Matlab，这里就给出用Matlab来进行面板数据的混合效应建模实例，特别是code解读

数据描述

该小组由美国48个州在1970年至1986年(17年)期间的年度观察结果组成。面板数据允许控制无法观察到的变量。例如，文化因素或企业间商业实践的差异; 或随时间而变化但不跨实体的变量(如国家政策、联邦法规、国际协议等)。也就是说，它解释了个体的异质性。该数据集由Munnell(1990)提供

1. 国内生产总值:国家经济分析局按州划分的国内生产总值

2. PUB_CAP:公共资本，包括(HWY, WATER, UTIL)

3. HWY: 高速公路和公路资本存量

4. WATER: 水务及污水设施资本存量

5. UTIL: 其他公共建筑、构筑物的资本存量

6. PVT_CAP: 根据美国国家经济分析局(Bureau of Economic Analysis)的私人资本存量估算

7. EMP: 非农就业人口，美国劳工统计局

8. UNEMP: 失业率，包括经济周期的影响，美国劳工统计局的所有美元数据都是数百万;就业数字以千为单位参考:

州和区域数据

1. GF = Gulf =AL, FL, LA, MS,

2. MW = Midwest =IL, IN, KY, Ml, MN, OH, Wl,

3. MA = Mid Atlantic =DE, MD, NJ, NY, PA, VA,

4. MT = Mountain =CO, ID, MT, ND, SD, WY,

5. NE = New England=CT, ME, MA, NH, Rl, VT,

6. SO = South =GA, NC, SC, TN, WV, R,

7. SW = Southwest =AZ, NV, NM, TX, UT,

8. CN = Central =AK, IA, KS, MO, NE, OK,

9. WC = West Coast =CA, OR, WA.

和中国的区域划分很相似，例如东，中，西，然后各包括那些省份等

调入数据

clear

clc

close all

load PanelData

预处理数据

将状态、年和区域转换为类别

publicdata.STATE =categorical(publicdata.STATE);

% Compute log(GDP)

publicdata.log_GDP = log(publicdata.GDP);

publicdata.GDP = [];

拟合线性模型，按区域可视化状态生产总值

数据收集于1970年至1986年之间的每个州。当拟合线性模型时，截距表示第0年的GDP。这可能会导致异常低的截距，由较大的坡度补偿。这导致了坡度估计值与相应的截距估计值之间的高度负相关。我们可以通过集中数据来消除这种相关性。

publicdata.YRcentered = publicdata.YR -mean(unique(publicdata.YR));

lm = fitlm(publicdata,'log_GDP ~ 1 +YRcentered');

可视化数据和拟合模型

figure,

gscatter(publicdata.YR,publicdata.log_GDP,publicdata.STATE,[],'.',20);

hold on

plot(publicdata.YR,predict(lm),'k-','LineWidth',2)

title('Simple linear modelfit','FontSize',15)

xlabel('Year','FontSize',15),ylabel('log(GDP)','FontSize',15)

从图中可以明显看出，拟合一个忽略状态(或区域)特定信息的简单线性回归模型并不足以完全指定log_GDP和YR之间的关系。这意味着我们忽略了组内任何违反线性回归假设的误差相关性，从而导致有偏差的标准误差。

figure,

xlabel('OLS Residuals')

boxplot(lm.Residuals.Raw,publicdata.STATE,'orientation','horizontal')

按状态划分的OLS残差箱线图显示，状态内残差要么全部为正，要么全部为负，且不以0为中心。这违反了OLS的基本假设，因此推断可能无效。

每组拟合单独的线性模型，并可视化参数的区间

figure,

plotIntervals(publicdata,'log_GDP ~ 1 +YRcentered','STATE')

单个置信区间清楚地表明，由于置信区间不重叠，因此需要一个随机效应来解释截距和斜率在州与州之间的变异性。

为每个状态拟合一个带有虚拟变量的线性模型

一种只使用固定效应来合并群体特定效应的方法是，引入状态的虚拟变量以及状态与年份之间的交互作用，如下所示

lm_group = fitlm(publicdata,'log_GDP ~ 1+ YRcentered + STATE + YRcentered:STATE');

% Display the number of predictors andnumber of estimated coefficients.

disp(' ')

disp(['Number of Variables:',num2str(lm_group.NumPredictors)])

disp(['Number of Estimated Coefficients:',...

   num2str(lm_group.NumEstimatedCoefficients)])

Number of Variables:2

Number of Estimated Coefficients:96

从图中可以看出，这个模型比之前的模型要好得多。只有当状态中有较少的能级和大量的观察时，它才是理想的。这种方法有几个缺点。尽管固定效应模型解释了状态效应，但它并没有提供有用的数据表示，因为它只对数据中的特定样本建模，而没有对总体进行概括。其次，模型中自由参数的数量随状态能级的数量线性增加，导致自由度的损失。在处理较小的观测集时，这可能是个问题。

拟合随机截距模型

让我们探索一个混合效应模型，在这个模型中，我们允许每个组的截距随机变化，在本例中是每个状态的截距随机变化。截距和YRcentered固定效应

clc

lme_intercept =fitlme(publicdata,'log_GDP ~ 1 + YRcentered + (1|STATE)');

状态随机效应项的标准差不包括0(下:0.82594，上:1.2324)，说明随机效应项具有显著性。

拟合随机截距和斜率模型

clc

lme_intercept_slope =fitlme(publicdata,...

   'log_GDP ~ 1 + YRcentered + (1+YRcentered|STATE)');

[~,~,rEffects] =randomEffects(lme_intercept_slope);

figure,

scatter(rEffects.Estimate(1:2:end),rEffects.Estimate(2:2:end))

title('Random Effects','FontSize',15)

xlabel('Intercept','FontSize',15)

ylabel('Slope','FontSize',15)

% The estimated column in therandom-effects table shows the estimated best

% linear unbiased predictor (BLUP) vectorof random effects.

用似然比检验比较随机斜率和截距模型

似然比检验(LRT)统计量可用于确定较复杂的模型是否明显优于较简单的模型。在这里，我们将随机截距模型与随机斜率和截距模型进行比较，以确定引入随机效应对边坡的意义。

compare(lme_intercept,lme_intercept_slope, 'CheckNesting',true)

ans =

    THEORETICAL LIKELIHOODRATIO TEST

    Model                  DF    AIC       BIC        LogLik    LRStat

    lme_intercept          4    -1550.5    -1531.7    779.24         

    lme_intercept_slope    6    -2085.4    -2057.2    1048.7   538.94

    deltaDF    pValue

    2          0    

小的p值表明，第二种模型解释了随机效应之间可能存在的相关性，明显优于第一种模型。“CheckNesting”标志尝试检查第一个模型是否嵌套在第二个模型中。

如果你不能简单地用公式描述你的模型，你可以创建设计矩阵来定义固定和随机的效果，然后调用fitlmematrix;

y = X*beta + Z*b + e

对于具有随机截距和斜率的模型，我们可以重现出结果，并且它们之间可能存在相关性

% Response (Dependent Variable)

y = publicdata.log_GDP;

% Fixed Effects Design Matrix

X = [ones(height(publicdata),1),publicdata.YRcentered];

% Random Effects Design Matrix

Z = [ones(height(publicdata),1),publicdata.YRcentered];

% Grouping Variable

G = publicdata.STATE;

% The random-effects design matrix canoptionally include the dummy coded

% variable for the group. If howevergroup 'G' is provided, fitlmematrix

% automatically does this for you.

% Fit the model

lme_matrix =fitlmematrix(X,y,Z,G,'CovariancePattern','Diagonal');

预测州GDP

states ={'IL','NJ','MA','CT','UT','NV','VT','WY'};

unbalancedStates = {'IL','CT','WY'};

missingYears = 1970:1983;

compareForecasts(publicdata, states,unbalancedStates, missingYears)

图中显示了两种方法对GDP预测的比较。实心圆表示用于拟合模型的观测值，空圆表示从模型中排除的观测值，但表示用于对预测进行可视化验证的观测值。第一个图显示了一个固定效应模型的预测，该模型使用州和州:年的虚拟变量来合并组特定的效应。第二幅图展示了随机斜率和截距模型对LME模型的预测。预测的置信区间表明，即使存在缺失的数据，LME模型也能提供更有信心的预测。对WY的进一步研究表明，由于可用的观测数据较少，固定效应模型显示，随着时间的推移，GDP预测会下降。另一方面，LME模型从整体上把握了数据的总体趋势，并预测了GDP的增长，这与缺失数据所显示的长期趋势相吻合。

最后是我自己练习的一个主程序代码

% Fixed Effects Panel Model withConcurrent Correlation

% Load sample data.

cd(matlabroot)

cd('help/toolbox/stats/examples')

load('panelData')

% Define variables

y = panelData.Growth;

city = panelData.City;

year = panelData.Year;

x = panelData.Employ;

% Plot data grouped by category.

figure(1)

set(0,'defaultfigurecolor','w')

boxplot(y,city)

set(gca,'FontName','Times NewRoman','FontSize',18);

xlabel('City','Fontsize',24)

saveas(gcf, 'zcpanle01.jpg')

% Plot data grouped by a differentcategory.

figure(2)

set(0,'defaultfigurecolor','w')

boxplot(y,year)

set(gca,'FontName','Times NewRoman','FontSize',18);

xlabel('Year','Fontsize',24)

saveas(gcf, 'zcpanle02.jpg')

% Format response data.

n = 6; d = 8;

Y = reshape(y,n,d);

% Format design matrices.

K = 7; N = n*d;

X = cell(n,1);

for i = 1:n

   x0 = zeros(d,K-1);

   x0(:,i) = 1;

   X{i} = [x0,x(i:n:N)];

end

% Fit the model.

[b,sig,E,V] =mvregress(X,Y,'algorithm','cwls');

% Plot fitted model.

xx = linspace(min(x),max(x));

axx = repmat(b(1:K-1),1,length(xx));

bxx = repmat(b(K)*xx,n,1);

yhat = axx + bxx;

figure(3)

set(0,'defaultfigurecolor','w')

hPoints = gscatter(x,y,year);

hold on

hLines = plot(xx,yhat);

set(gca,'FontName','Times NewRoman','FontSize',18);

for i=1:n 

  set(hLines(i),'color',get(hPoints(i),'color'));

end

hold off

saveas(gcf, 'zcpanle03.jpg')

% Residual correlation.

figure(4)

set(0,'defaultfigurecolor','w')

gscatter(year,E(:),city)

set(gca,'FontName','Times NewRoman','FontSize',18);

ylabel('Residuals','Fontsize',24)

saveas(gcf, 'zcpanle04.jpg')

% Panel corrected standard errors.

XX = cell2mat(X);

S = kron(eye(n),sig);

Vpcse = inv(XX'*XX)*XX'*S*XX*inv(XX'*XX);

se = sqrt(diag(Vpcse))