自习时偶然发现的一段代码 需要的自取哈—帮原作者造福作业小白/寻找原作者/如果有侵权行为立马告知、速删 #经济学代码#MATLAB#代码自取#模型预测代码

第一题：

1. S1是一个季节性模型，包括第四季度的虚拟变量和线性趋势；S2是一个季节性模型，允许所有季度和线性趋势；S3是一个月度模型，包括第12个月的虚拟变量和线性趋势；S4是一个月度模型，允许所有月份和线性趋势。

（2）

% import Australia employment and GDP data

Data = importdata("AUSEmploy2022.csv");

data = Data.data;

textdata = Data.textdata;

% wash the data and generate the dummy variable of month and quarter

dt = datetime(textdata(2:end,1),'InputFormat','dd/MM/uuuu');

Year = year(dt);

Quarter = quarter(dt);

Month = month(dt);

y = data(:,1);

x = data(:,2);

T = length(y);

M = zeros(T,12);

D = zeros(T,4);

for m = 1:12

    M(:,m) = (Month==m);

end % month dummy

for q = 1:4

    D(:,q) = (Quarter==q);

end % quarter dummy

% S1

X1 = [ones(T,1) (1:T)' D(:,4)];

beta1 = (X1'*X1)\(X1'*y);

yhat1 = X1*beta1;

MSE1 = mean((y-yhat1).^2);

AIC1 = T*MSE1 + 3*2;

BIC1 = T*MSE1 + 3*log(T);

% S2

X2 = [(1:T)' D];

beta2 = (X2'*X2)\(X2'*y);

yhat2 = X2*beta2;

MSE2 = mean((y-yhat2).^2);

AIC2 = T*MSE2 + 5*2;

BIC2 = T*MSE2 + 5*log(T);

% S3

X3 = [ones(T,1) (1:T)' M(:,12)];

beta3 = (X3'*X3)\(X3'*y);

yhat3 = X3*beta3;

MSE3 = mean((y-yhat3).^2);

AIC3 = T*MSE3 + 3*2;

BIC3 = T*MSE3 + 3*log(T);

% S4

X4 = [(1:T)' M];

beta4 = (X4'*X4)\(X4'*y);

yhat4 = X4*beta4;

MSE4 = mean((y-yhat4).^2);

AIC4 = T*MSE4 + 13*2;

BIC4 = T*MSE4 + 13*log(T);

Model = {'S1';'S2';'S3';'S4'};

MSE = [MSE1;MSE2;MSE3;MSE4];

AIC = [AIC1;AIC2;AIC3;AIC4];

BIC = [BIC1;BIC2;BIC3;BIC4];

table(Model,MSE,AIC,BIC) % tabulate the MSE AIC BIC of each model

ans = 4×4 table

	Model	MSE	AIC	BIC
1	'S1'	2.2717e+10	7.3829e+12	7.3829e+12
2	'S2'	2.2589e+10	7.3414e+12	7.3414e+12
3	'S3'	2.1896e+10	7.1163e+12	7.1163e+12
4	'S4'	2.0321e+10	6.6044e+12	6.6044e+12

beta1' %  parameters of S1

ans = 1×3

106 ×

    7.8529    0.0163    0.0575

beta2' %  parameters of S2

ans = 1×5

106 ×

    0.0163    7.8364    7.8681    7.8549    7.9106

beta3' %  parameters of S3

ans = 1×3

106 ×

    7.8559    0.0163    0.1375

beta4' %  parameters of S4

ans = 1×13

106 ×

    0.0163    7.7562    7.8679    7.8870    7.8665    7.8650    7.8720

    7.8730    7.7902    7.9004    7.8599    7.8778    7.9931

（3）

T0 = 50;

h1 = 1; h2 = 5;

ytph1 = y(T0+h1:end); ytph2 = y(T0+h2:end);

syhat1 = zeros(T-T0-h1+1,4); syhat2 = zeros(T-T0-h2+1,4);

% 1-step-ahead forcast

for t = T0:T-h1

    XX1 = X1(1:t,:);

    XX2 = X2(1:t,:);

    XX3 = X3(1:t,:);

    XX4 = X4(1:t,:);

    yy = y(1:t);

    beta1 = (XX1'*XX1)\(XX1'*yy);

    beta2 = (XX2'*XX2)\(XX2'*yy);

    beta3 = (XX3'*XX3)\(XX3'*yy);

    beta4 = (XX4'*XX4)\(XX4'*yy);

    syhat1(t-T0+1,1) = X1(t+h1,:)*beta1;

    syhat1(t-T0+1,2) = X2(t+h1,:)*beta2;

    syhat1(t-T0+1,3) = X3(t+h1,:)*beta3;

    syhat1(t-T0+1,4) = X4(t+h1,:)*beta4;

end

MSFE1 = mean((syhat1-ytph1).^2); MSFE1_Q1 = MSFE1;

% 5-step-ahead forcast

for t = T0:T-h2

    XX1 = X1(1:t,:);

    XX2 = X2(1:t,:);

    XX3 = X3(1:t,:);

    XX4 = X4(1:t,:);

    yy = y(1:t);

    beta1 = (XX1'*XX1)\(XX1'*yy);

    beta2 = (XX2'*XX2)\(XX2'*yy);

    beta3 = (XX3'*XX3)\(XX3'*yy);

    beta4 = (XX4'*XX4)\(XX4'*yy);

    syhat2(t-T0+1,1) = X1(t+h2,:)*beta1;

    syhat2(t-T0+1,2) = X2(t+h2,:)*beta2;

    syhat2(t-T0+1,3) = X3(t+h2,:)*beta3;

    syhat2(t-T0+1,4) = X4(t+h2,:)*beta4;

end

MSFE2 = mean((syhat2-ytph2).^2); MSFE2_Q1 = MSFE2;

Table = table(Model,MSFE1',MSFE2');

Table.Properties.VariableNames = {'Model','MSFE(1-step)','MSFE(5-step)'}

Table = 4×3 table

	Model	MSFE(1-step)	MSFE(5-step)
1	'S1'	2.6831e+10	3.0950e+10
2	'S2'	2.3225e+10	2.3456e+10
3	'S3'	2.5898e+10	2.9920e+10
4	'S4'	2.5816e+10	2.9224e+10

（4）

figure

subplot(2,1,1)

plot(dt(T0+h1:end),[ytph1 syhat1]),datetick('x','yyyy'),legend({'Origin','S1','S2','S3','S4'},'Location','eastoutside'),title('1-step-ahead forecast')

subplot(2,1,2)

plot(dt(T0+h2:end),[ytph2 syhat2]),datetick('x','yyyy'),legend({'Origin','S1','S2','S3','S4'},'Location','eastoutside'),title('5-step-ahead forecast')

根据问题b的结果，S4具有更好的拟合优度，但是根据问题c的结果，S2是样本外预测的最佳模型。因此，S4是过度拟合的，而S2是最合适的模型。

 

第二题

（1）

k = [3 5 3 13]'; % the number of parameters of each model

sigma2 = MSE*T./k;

table(Model,sigma2)

ans = 4×2 table

	Model	sigma2
1	'S1'	2.4610e+12
2	'S2'	1.4683e+12
3	'S3'	2.3721e+12
4	'S4'	5.0803e+11

（2）

idx1 = find(dt(T0+h1:end)==datetime(2021,6,1));

lower1 = syhat1(idx1,:) - sqrt(sigma2')*norminv(0.975);

upper1 = syhat1(idx1,:) + sqrt(sigma2')*norminv(0.975);

idx2 = find(dt(T0+h2:end)==datetime(2021,6,1));

lower2 = syhat2(idx2,:) - sqrt(sigma2')*norminv(0.975);

upper2 = syhat2(idx2,:) + sqrt(sigma2')*norminv(0.975);

Table = table(Model,[lower1' upper1'],[lower2' upper2']);

Table.Properties.VariableNames = {'Model','95% interval(1-step)','95% interval(5-step)'}

Table = 4×3 table

	Model	95% interval(1-step)		⋯
1	'S1'	9.9565e+06	1.6106e+07
2	'S2'	1.0671e+07	1.5421e+07
3	'S3'	1.0016e+07	1.6053e+07
4	'S4'	1.1647e+07	1.4441e+07

（3）

X1 = [ones(T,1) (1:T)' D(:,4)];

beta1 = (X1'*X1)\(X1'*y);

X2 = [(1:T)' D];

beta2 = (X2'*X2)\(X2'*y);

X3 = [ones(T,1) (1:T)' M(:,12)];

beta3 = (X3'*X3)\(X3'*y);

X4 = [(1:T)' M];

beta4 = (X4'*X4)\(X4'*y);

% in 2022-06

XX1 = [1 T+5 0];

XX2 = [T+5 0 1 0 0];

XX3 = [1 T+5 0];

XX4 = [T+5 0 0 0 0 0 1 0 0 0 0 0 0];

ypred1 = XX1*beta1;

ypred2 = XX2*beta2;

ypred3 = XX3*beta3;

ypred4 = XX4*beta4;

ypred = [ypred1 ypred2 ypred3 ypred4]';

Prob = 1 - normcdf((13.5e6-ypred)./sqrt(sigma2));

table(Model,Prob)

Table = 4×2 table

	Model	Prob
1	'S1'	0.4307
2	'S2'	0.4154
3	'S3'	0.4302
4	'S4'	0.3606

（4）

使用20个滞后的Ljung-Box Q-test来测试每个模型的误差的自相关性。

error1 = y - X1*beta1;

h1 = lbqtest(error1);

error2 = y - X2*beta2;

h2 = lbqtest(error2);

error3 = y - X3*beta3;

h3 = lbqtest(error3);

error4 = y - X4*beta4;

h4 = lbqtest(error4);

H = [h1 h2 h3 h4]';

table(Model,H)

ans = 4×2 table

	Model	H
1	'S1'	1
2	'S2'	1
3	'S3'	1
4	'S4'	1

使用残差自相关的Ljung-Box Q-test，四个模型的误差都拒绝H0:一系列残差没有表现出自相关。所以et是iid对于这个数据序列是一个不合理的假设。

第三题

T0 = 50;

h1 = 1; h2 = 5;

ytph1 = y(T0+h1:end); ytph2 = y(T0+h2:end);

syhat1 = zeros(T-T0-h1+1,4); syhat2 = zeros(T-T0-h2+1,4);

% 1-step-ahead forcast

for t = T0:T-h1

    XXX1 = [X1(1+h1:t,:) x(1:t-h1)];

    XXX2 = [X2(1+h1:t,:) x(1:t-h1)];

    XXX3 = [X3(1+h1:t,:) x(1:t-h1)];

    XXX4 = [X4(1+h1:t,:) x(1:t-h1)];

    yyy = y(1+h1:t);

    beta1 = (XXX1'*XXX1)\(XXX1'*yyy);

    beta2 = (XXX2'*XXX2)\(XXX2'*yyy);

    beta3 = (XXX3'*XXX3)\(XXX3'*yyy);

    beta4 = (XXX4'*XXX4)\(XXX4'*yyy);

    syhat1(t-T0+1,1) = [X1(t+h1,:) x(t)]*beta1;

    syhat1(t-T0+1,2) = [X2(t+h1,:) x(t)]*beta2;

    syhat1(t-T0+1,3) = [X3(t+h1,:) x(t)]*beta3;

    syhat1(t-T0+1,4) = [X4(t+h1,:) x(t)]*beta4;

end

MSFE1 = mean((syhat1-ytph1).^2); MSFE1_Q3 = MSFE1;

% 5-step-ahead forcast

for t = T0:T-h2

    XXX1 = [X1(1+h2:t,:) x(1:t-h2)];

    XXX2 = [X2(1+h2:t,:) x(1:t-h2)];

    XXX3 = [X3(1+h2:t,:) x(1:t-h2)];

    XXX4 = [X4(1+h2:t,:) x(1:t-h2)];

    yyy = y(1+h2:t);

    beta1 = (XXX1'*XXX1)\(XXX1'*yyy);

    beta2 = (XXX2'*XXX2)\(XXX2'*yyy);

    beta3 = (XXX3'*XXX3)\(XXX3'*yyy);

    beta4 = (XXX4'*XXX4)\(XXX4'*yyy);

    syhat2(t-T0+1,1) = [X1(t+h2,:) x(t)]*beta1;

    syhat2(t-T0+1,2) = [X2(t+h2,:) x(t)]*beta2;

    syhat2(t-T0+1,3) = [X3(t+h2,:) x(t)]*beta3;

    syhat2(t-T0+1,4) = [X4(t+h2,:) x(t)]*beta4;

end

MSFE2 = mean((syhat2-ytph2).^2); MSFE2_Q3 = MSFE2;

Model = {'S1 prime','S2 prime','S3 prime','S4 prime'}';

Table = table(Model,MSFE1',MSFE2');

Table.Properties.VariableNames = {'Model','MSFE(1-step)','MSFE(5-step)'}

Table = 4×3 table

	Model	MSFE(1-step)	MSFE(5-step)
1	'S1 prime'	2.3373e+10	3.0711e+10
2	'S2 prime'	2.3440e+10	2.3772e+10
3	'S3 prime'	2.2442e+10	2.9756e+10
4	'S4 prime'	2.1814e+10	2.9009e+10

figure

subplot(2,1,1)

plot(dt(T0+h1:end),[ytph1 syhat1]),datetick('x','yyyy'),legend({'Origin','S1 prime','S2 prime','S3 prime','S4 prime'},'Location','eastoutside'),title('1-step-ahead forecast')

subplot(2,1,2)

plot(dt(T0+h2:end),[ytph2 syhat2]),datetick('x','yyyy'),legend({'Origin','S1 prime','S2 prime','S3 prime','S4 prime'},'Location','eastoutside'),title('5-step-ahead forecast')

在一步预测中，S4 为主是最好的模型，其最小MFSE为218.14亿。在5步预测中，S2 为主是最好的模型，其最小MFSE为237.72亿。

 

我们可以在5步预测中得出结论，S4素数是过度拟合的。

第四题

（1）

alpha = 0.3; beta = 0.3; gamma = 0.3;

T0 = 50;

h1 = 1; h2 = 5;

s = 4;

syhat1 = zeros(T-h1-T0+1,1); syhat2 = zeros(T-h2-T0+1,1);

ytph1 = y(T0+h1:end); ytph2 = y(T0+h2:end);

St1 = zeros(T-min(h1,s),1); St2 = zeros(T-min(h2,s),1);

Lt = mean(y(1:s)); bt = 0; St1(1:s) = y(1:s) - Lt;

for t = s+1:T-h1

    newLt = alpha*(y(t)-St1(t-s)) + (1-alpha)*(Lt+bt);

    newbt = beta*(newLt-Lt) + (1-beta)*bt;

    St1(t) = gamma*(y(t)-newLt) + (1-gamma)*St1(t-s);

    yhat = newLt + h1*newbt + St1(t+h1-s);

    Lt = newLt; bt = newbt;

    if t >= T0

        syhat1(t-T0+1,:) = yhat;

    end

end

MSFE1_HW = mean((ytph1-syhat1).^2);

Lt = mean(y(1:s)); bt = 0; St2(1:s) = y(1:s) - Lt;

for t = s+1:T-h2

    newLt = alpha*(y(t)-St2(t-s)) + (1-alpha)*(Lt+bt);

    newbt = beta*(newLt-Lt) + (1-beta)*bt;

    St2(t) = gamma*(y(t)-newLt) + (1-gamma)*St2(t-s);

    yhat = newLt + h2*newbt + St2(t+h2-s);

    Lt = newLt; bt = newbt;

    if t >= T0

        syhat2(t-T0+1,:) = yhat;

    end

end

MSFE2_HW = mean((ytph2-syhat2).^2);

Table = table(MSFE1_HW,MSFE2_HW);

Table.Properties.VariableNames = {'MSFE(1-step)','MSFE(5-step)'}

Table = 1×2 table

	MSFE(1-step)	MSFE(5-step)
1	1.4262e+10	3.4288e+10

（2）

params0 = [0.3 0.3 0.3];

[bestparams1,bestMSFE1_HW] = fmincon(@(params) HWMFSE(y,h1,s,T0,params),params0,[],[],[],[],[0 0 0],[1 1 1]);

Local minimum possible. Constraints satisfied.

fmincon stopped because the size of the current step is less than

the value of the step size tolerance and constraints are

satisfied to within the value of the constraint tolerance.

<stopping criteria details>

[bestparams2,bestMSFE2_HW] = fmincon(@(params) HWMFSE(y,h2,s,T0,params),params0,[],[],[],[],[0 0 0],[1 1 1]);

Local minimum possible. Constraints satisfied.

fmincon stopped because the size of the current step is less than

the value of the step size tolerance and constraints are

satisfied to within the value of the constraint tolerance.

<stopping criteria details>

reduction1 = (MSFE1 -bestMSFE1_HW)/MSFE1*100; reduction2 = (MSFE2-bestMSFE2_HW)/MSFE2*100;

rownames = {'alpha','beta','gamma','MFSE','reduction(%)'}';

variablenames = {' ','best 1-step-ahead forecast','best 5-step-ahead forecast'};

Table = table(rownames,[bestparams1 bestMSFE1_HW reduction1]',[bestparams2 bestMSFE2_HW reduction2]');

Table.Properties.VariableNames = variablenames

Table = 5×3 table

		best 1-step-ahead forecast	best 5-step-ahead forecast
1	'alpha'	0.4734	0.3546
2	'beta'	0.0143	0.0165
3	'gamma'	0.1284	0.0023
4	'MFSE'	1.0911e+10	2.0629e+10
5	'reduction(%)'	52.1197	27.7900

与基线相比，最佳硬件模型的MFSEs在最佳1步预测中减少了52.1197%，在最佳5步预测中减少了27.7900%。

第五题

rownames = {'S1','S2','S3','S4','S1 prime','S2 prime','S3 prime','S4 prime','HW','bestHW'}';

MSFE1 = [MSFE1_Q1 MSFE1_Q3 MSFE1_HW bestMSFE1_HW]';

MSFE2 = [MSFE2_Q1 MSFE2_Q3 MSFE2_HW bestMSFE2_HW]';

figure

subplot(2,1,1)

plot(MSFE1),xticklabels(rownames),ylabel('MSFE'),title('1-step-ahead forecast')

subplot(2,1,2)

plot(MSFE2),xticklabels(rownames),ylabel('MSFE'),title('5-step-ahead forecast')

显然，最好的HW是我们用过的所有型号中最好的型号。具体来说，S4主要在1步提前预测和S2在5步提前预测表现良好，没有HW模型。

 

第六题

1. 和（2）

我决定使用SARIMAX模型和外生变量澳大利亚月度GDP。

根据（2 ）d中的分析，S1-S4模型的误差表现出自相关性。所以我选择SARIMAX模型来捕捉自相关性、季节性和外生因素的特征。

（3）

模型的设置如下:

 

 

1. 和（5）

Md = arima('ARLags',1,'D',1,'MALags',1,'Seasonality',12);

idxpre = 1:Md.P;

syhat1 = zeros(T-T0-h1+1,1); syhat2 = zeros(T-T0-h2+1,1);

ytph1 = y(T0+h1:end); ytph2 = y(T0+h2:end);

for t = T0:T-h1

    idxest = Md.P+1:t;

    EstMd = estimate(Md,y(idxest),'Y0',y(idxpre),'X',x(idxest),'Display','off');

    yf = forecast(EstMd,h1,y(idxest(end-Md.P+1:end)),'XF',x(t+1:t+h1));

    syhat1(t-T0+1) = yf(end);

end

myMSFE1 = mean((syhat1-ytph1).^2);

for t = T0:T-h2

    idxest = Md.P+1:t;

    EstMd = estimate(Md,y(idxest),'Y0',y(idxpre),'X',x(idxest),'Display','off');

    yf = forecast(EstMd,h2,y(idxest(end-Md.P+1:end)),'XF',x(t+1:t+h2));

    syhat2(t-T0+1) = yf(end);

end

myMSFE2 = mean((syhat2-ytph2).^2);

Table = table(myMSFE1,myMSFE2);

Table.Properties.VariableNames = {'MSFE(1-step)','MSFE(5-step)'}

Table = 1×2 table

	MSFE(1-step)	MSFE(5-step)
1	8.9343e+09	4.9231e+10

rownames = {'S1','S2','S3','S4','S1 prime','S2 prime','S3 prime','S4 prime','HW','bestHW','SARIMAX'}';

MSFE1 = [MSFE1_Q1 MSFE1_Q3 MSFE1_HW bestMSFE1_HW myMSFE1]';

MSFE2 = [MSFE2_Q1 MSFE2_Q3 MSFE2_HW bestMSFE2_HW myMSFE2]';

figure

subplot(2,1,1)

plot(MSFE1),xticklabels(rownames),ylabel('MSFE'),title('1-step-ahead forecast')

subplot(2,1,2)

plot(MSFE2),xticklabels(rownames),ylabel('MSFE'),title('5-step-ahead forecast')

 

根据上图，SARIMAX模型在之前使用的所有模型中在1步预测中表现最好，但在5步预测中表现最差。

第七题

为了更合理地设置φ1和φ2的值，采用两步OLS法计算所有参数的先验。

XXXX = [ones(T-1,1) (2:T)' y(1:T-1)];

yyyy = y(2:T);

beta = (XXXX'*XXXX)\(XXXX'*yyyy);

mu = beta(1); a = beta(2); rho = beta(3);

e = y - [sum(beta(1:2)); XXXX*beta];

XXXX = [e(2:T-1) e(1:T-2)];

yyyy = e(3:T);

beta = (XXXX'*XXXX)\(XXXX'*yyyy);

phi = beta;

u = yyyy - XXXX*beta;

sigma2 = mean(u'*u);

    Use MLE for

 

format long

params0 = [mu a rho sigma2];

params = fmincon(@(params) ARX1LL(y,phi,params),params0);

Problem appears unbounded.

fmincon stopped because the objective function value is less than

the value of the objective function limit and constraints

are satisfied to within the value of the constraint tolerance.

<stopping criteria details>

params = 1×4

1012 ×

   0.000002204395936  -0.000001588627601  -0.096718249297454 ⋯

mu = params(1); a = params(2); rho = params(3); sigma2 = params(4);

table(mu,a,rho,sigma2)

ans = 1×4 table

	mu	a	rho	sigma2
1	2.2044e+06	-1.5886e+06	-9.6718e+10	3.4294e+12

作业中使用的两个公式：

function out = HWMFSE(y,h,s,T0,params)

T = length(y);

alpha = params(1); beta = params(2); gamma = params(3);

syhat = zeros(T-h-T0+1,1);

ytph = y(T0+h:end);

St = zeros(T-min(h,s),1);

Lt = mean(y(1:s)); bt = 0; St(1:s) = y(1:s) - Lt;

for t = s+1:T-h

    newLt = alpha*(y(t)-St(t-s)) + (1-alpha)*(Lt+bt);

    newbt = beta*(newLt-Lt) + (1-beta)*bt;

    St(t) = gamma*(y(t)-newLt) + (1-gamma)*St(t-s);

    yhat = newLt + h*newbt + St(t+h-s);

    Lt = newLt; bt = newbt;

    if t >= T0

        syhat(t-T0+1,:) = yhat;

    end

end

out = mean((ytph-syhat).^2);

end

function out = ARX1LL(y,phi,params)

T = length(y);

phi1 = phi(1); phi2 = phi(2);

mu = params(1); a = params(2); rho = params(3); sigma2 = params(4);

L = zeros(T);

for i = 2:T

    L(i,i-1) = 1;

end

I = eye(T);

e = (I-rho*L)*y - mu - a*(1:T)';

u = (I - phi1*L - phi2*L^2)*e;

out = -pi/2*log(2*pi*sigma2) - 1/(2*sigma2)*(u'*u);

end