clear all;

% 工作路径
userpath('C:\Users\Surface\Desktop'); 
% 读取数据
[data,name] = xlsread('data',1); % 原数据为原始面板数据，需要转化为空间面板数据

% 读取几个权重矩阵
weight1 = xlsread('华东各省会城市距离',1);            % 基于距离关系的空间权重矩阵
weight2 = xlsread('华东省区公路邻域',1);              % 基于邻域关系的空间权重矩阵
weight3t = xlsread('华东2014GDP计算的经济距离',1);    % 基于GDP的经济距离关系的空间权重矩阵
weight3 = fanzhuan(weight3t);                        % 将上三角矩阵反转为实对称矩阵，需要fanzhuan包
% 空间权重标准化,函数需要fanzhuan包
W1 = hbzh(weight1);
W2 = lyhbzh(weight2);
W3 = hbzh(weight3);
W = 1/3*W1+1/3*W2+1/3*W3    % 加权构建复合矩阵

% 问题的维度
T = 40; % 时期数量 
N = 7; % 地区数量

% 调整数据结构,从原始数据转化为空间面板数据
temp =  zeros(T,7);  % 构建存储面板数据的0矩阵
for i = 1:T         % 第i年至第40年
    
    shh = data(i,[2:8]);  % SHH上海当年的全部数据
    shd = data(i,[9:15]); % SHD山东当年的全部数据
    ah = data(i,[16:22]);  % AH 安徽当年的全部数据
    js = data(i,[23:29]);  % JS 江苏当年的全部数据
    zhj = data(i,[30:36]); % ZHJ浙江当年的全部数据
    jx = data(i,[37:43]);  % JX 江西当年的全部数据
    fj = data(i,[44:50]);  % FJ 福建当年的全部数据
    temp{(7*i-6):(7*i),:) = vertcat(shh,shd,ah,js,zhj,jx,fj);  % 将每一年的数据贴进准备好的0矩阵
    
    if i>40
        break
    end
end

% 0 OLS模型
% 提取变量
y = temp(:,[1]);           % 华东各地区GDP
x_ols = temp(:,[2:7]); 
% 变量名
vnames=strvcat('gdp','renkou','gsh','gzh','tz','xf','jk'); 
% 因变量：GDP（亿元），自变量：人口(万人)，财政收入（亿元），财政支出（亿元），投资（亿元），消费（亿元），进出口（亿元）


% 建模
results_ols = ols(y,x_ols); % 进行最小二乘法回归
prt_reg(results_ols,vnames,1); 

I = eye(T);W1 = kron(I,W);
res = moran(y,x_ols,W1); % moran 检验
prt(res); %moran值较大，在0.01显著性水平下，依据空间分布的特征显著

% ---------------------------------------------------------------------
% 1 SAR模型
% 提取变量：2017年数据
y = temp(274:280,[1]);         
x_sar = temp(274:280,[2,3,4,5,6,7]);      
% 因变量：GDP（亿元），自变量：人口(万人)，财政收入（亿元），财政支出（亿元），投资（亿元），消费（亿元），进出口（亿元）

% 建模
results_sar = sar(y,x_sar,W);
prt(results_sar,vnames);

results_sar.tstat        % t检验统计量的值
results_sar.rho          % 空间自回归系数rho
results_sar.beta         % 模型估计的beta值


% ---------------------------------------------------------------------
% 2 SEM模型
% y = beta*x + mu ; mu = rho*W*mu + e 
% 选择变量：2017年数据
y = temp(274:280,[1]); 
x_sem = temp(274:280,[2:7]); 

results_sem = sem(y,x_sem,W);
prt(results_sem,vnames);

results_sem.tstat    % 渐进t检验统计量的值(最后输入为rho空间自回归系数)
results_sem.rho      % (p above)
results_sem.beta     % 模型估计的beta值


% ---------------------------------------------------------------------
% 3 SDM空间杜宾模型
% y = XB + u (optional) + v (optional) + s,  s = p*W*s + e
% 选择变量:2017年数据
y = temp(274:280,[1]);  
x_sdm = temp(274:280,[2:4]); % 空间相关性影响变量
vnames1=strvcat('gdp','renkou','gsh','gzh'); 

% 建模
results_sdm = sdm(y,x_sdm,W); 
prt(results_sdm,vnames1);

results_sdm.tstat        % t检验统计量的值
results_sdm.rho          % 最后输入为rho空间自回归系数
results_sdm.beta         % 模型估计的beta值

% 4 空间面板数据模型

% 4.1 空间相关性检验
%LM检验，大于6.635拒绝原假设（不存在空间自相关），即存在误差项空间自相关
LMsarsem_panel(results_sar,W,y,x_sar) 
% 结果表明0.1显著性水平下，空间固定效应LR检验通过，空间固定效应的空间滞后的LM检验通过
% 时间固定效应的的空间误差的LM检验不通过,时间固定效应的LR检验不通过
% 故选择时间效应固定模型info.model=2
info.model=2;

% 4.2 提取变量：30年数据
y = temp(:,[1]);         
x_sar = temp(:,[2,3,4,5,6,7]);      
% 因变量：GDP（亿元），自变量：人口(万人)，财政收入（亿元），财政支出（亿元），投资（亿元），消费（亿元），进出口（亿元）

% 4.3 建模
results_general = sar_panel_FE(y,x_sar,W,T,info); 
prt_spnew(results_general,vnames);

results_general.corr2     % 模型的拟合优度调整的r方
results_general.tstat     % t检验统计量的值
results_general.rho       % 最后输入为rho空间自回归系数
results_general.beta      % 模型估计的beta值