%清理运行环境
clear;
clc;
close all;
addpath('./psopt');
global itercount;
itercount = 0;

%读取沪深300指数2005年-2016年上半年的历史数据
data = xlsread('data.xls', 'Sheet1', 'B1:B365');
N = numel(data);
y = data(4 : end);
x = [data(3 : N - 1), data(2 : N - 2), data(1 : N - 3)];

%构造随机微分方程
func = @(xx) myfunc(xx, [y, x]);
%估计随机微分方程模型的参数
options = psooptimset('ConstrBoundary', 'reflect', 'DemoMode', 'pretty', 'HybridFcn', @fminunc, 'TolCon', 1e-6);
[xOpt,fval,exitflag,output,population,scores] = pso(func, 8, [], [], [], [], [], [], [], options);

%套算预测结果
x1 = xOpt(1);
x2 = xOpt(2);
x3 = xOpt(3);
x4 = xOpt(4);
y1 = xOpt(5);
y2 = xOpt(6);
y3 = xOpt(7);
y4 = xOpt(8);

xx = x;
I = ones(numel(y), 1);
xx1 = [I, xx];

%按照多项式规则构造随机微分方程的均值与标准差项
fx = xx1 * [x1; x2; x3; x4];
gx = xx1 * [y1; y2; y3; y4];

%使用10000次模拟随机微分方程的计算过程
rand('seed', 0);
numberOfRandom = 10000;
randomNumbers = normrnd(0, 1, numel(y), numberOfRandom);
tmp = repmat(gx, 1, numberOfRandom) .* randomNumbers;
tmp1 = sum(tmp, 2) / numberOfRandom;
py = fx + tmp1;
residuals = y - py;
RMSE = sqrt(mean((residuals ./ y) .^ 2));

rmpath('./psopt');

save result.mat;
system('shutdown -s');