几何布朗运动模拟 MATLAB实现

time = data.date;%抽取年份
close = data.close; %抽取收盘价
pos = time.Year == 2019;%选择2019年打标签
close_2019 = close(pos);%选择2019年收盘价

[ret,vol] = RetVolCaculation(data);
simup = MCSimulation(ret,vol,time,close);
resultplot(close_2019,simup);

function [ret,vol] = RetVolCaculation(data)
ret = zeros(14,1);
vol = zeros(14,1);
time_2005 = table2array(data(timerange...
    ('2005-01-04','2005-12-30'),'close'));
time_2006 = table2array(data(timerange...
    ('2006-01-04','2006-12-29'),'close'));
time_2007 = table2array(data(timerange...
    ('2007-01-04','2007-12-28'),'close'));
time_2008 = table2array(data(timerange...
    ('2008-01-04','2008-12-29'),'close'));
time_2009 = table2array(data(timerange...
    ('2009-01-05','2009-12-11'),'close'));
time_2010 = table2array(data(timerange...
    ('2010-01-04','2010-12-30'),'close'));
time_2011 = table2array(data(timerange...
    ('2011-01-04','2011-12-30'),'close'));
time_2012 = table2array(data(timerange...
    ('2012-01-04','2012-12-28'),'close'));
time_2013 = table2array(data(timerange...
    ('2013-01-04','2013-12-30'),'close'));
time_2014 = table2array(data(timerange...
    ('2014-01-02','2014-12-30'),'close'));
time_2015 = table2array(data(timerange...
    ('2015-01-05','2015-12-30'),'close'));
time_2016 = table2array(data(timerange...
    ('2016-01-05','2016-12-29'),'close'));
time_2017 = table2array(data(timerange...
    ('2017-01-03','2017-12-29'),'close'));
time_2018 = table2array(data(timerange...
    ('2018-01-02','2018-12-28'),'close'));


ret(1) = log(table2array(data({'2005-01-04'},'close')))...
    - log(table2array(data({'2005-12-30'},'close')));
ret(2) = log(table2array(data({'2006-01-04'},'close')))...
    - log(table2array(data({'2006-12-29'},'close')));
ret(3) = log(table2array(data({'2007-01-04'},'close')))...
    - log(table2array(data({'2007-12-28'},'close')));
ret(4) = log(table2array(data({'2008-01-04'},'close')))...
    - log(table2array(data({'2008-12-29'},'close')));
ret(5) =log(table2array(data({'2009-01-05'},'close')))...
    - log(table2array(data({'2009-12-11'},'close')));
ret(6) = log(table2array(data({'2010-01-04'},'close')))...
    - log(table2array(data({'2010-12-30'},'close')));
ret(7) = log(table2array(data({'2011-01-04'},'close')))...
    - log(table2array(data({'2011-12-30'},'close')));
ret(8) = log(table2array(data({'2012-01-04'},'close')))...
    - log(table2array(data({'2012-12-28'},'close')));
ret(9) =  log(table2array(data({'2013-01-04'},'close')))...
    - log(table2array(data({'2013-12-30'},'close')));
ret(10) =  log(table2array(data({'2014-01-02'},'close')))...
    - log(table2array(data({'2014-12-30'},'close')));
ret(11) =  log(table2array(data({'2015-01-05'},'close')))...
    - log(table2array(data({'2015-12-30'},'close')));
ret(12) =  log(table2array(data({'2016-01-05'},'close')))...
    - log(table2array(data({'2016-12-29'},'close')));
ret(13) =  log(table2array(data({'2017-01-03'},'close')))...
    - log(table2array(data({'2017-12-29'},'close')));
ret(14) =  log(table2array(data({'2018-01-02'},'close')))...
    - log(table2array(data({'2018-12-28'},'close')));


m = length(time_2005);
vol(1) = 0;
for i = 1:m-1
    vol(1) = vol(1) + (log(time_2005(i+1)) - log(time_2005(i)))^2;
end
vol(1) = sqrt(vol(1));

m = length(time_2006);
vol(2) = 0;
for i = 1:m-1
    vol(2) = vol(2) + (log(time_2006(i+1)) - log(time_2006(i)))^2;
end
vol(2) = sqrt(vol(2));

m = length(time_2007);
vol(3) = 0;
for i = 1:m-1
    vol(3) = vol(3) + (log(time_2007(i+1)) - log(time_2007(i)))^2;
end
vol(3) = sqrt(vol(3));

m = length(time_2008);
vol(4) = 0;
for i = 1:m-1
    vol(4) = vol(4)+(log(time_2008(i+1)) - log(time_2008(i)))^2;
end
vol(4) = sqrt(vol(4));

m = length(time_2009);
vol(5) = 0;
for i = 1:m-1
    vol(5) = vol(5) + (log(time_2009(i+1)) - log(time_2009(i)))^2;
end
vol(5) = sqrt(vol(5));

m = length(time_2010);
vol(6) = 0;
for i = 1:m-1
    vol(6) = vol(6)+(log(time_2010(i+1)) - log(time_2010(i)))^2;
end
vol(6) = sqrt(vol(6));

m = length(time_2011);
vol(7) = 0;
for i = 1:m-1
    vol(7) =vol(7)+ (log(time_2011(i+1)) - log(time_2011(i)))^2;
end
vol(7) = sqrt(vol(7));

m = length(time_2012);
vol(8) = 0;
for i = 1:m-1
    vol(8) = vol(8) + (log(time_2012(i+1)) - log(time_2012(i)))^2;
end
vol(8) = sqrt(vol(8));

m = length(time_2013);
vol(9) = 0;
for i = 1:m-1
    vol(9) = vol(9) +  (log(time_2013(i+1)) - log(time_2013(i)))^2;
end
vol(9) = sqrt(vol(9));

m = length(time_2014);
vol(10) = 0;
for i = 1:m-1
    vol(10) = vol(10) + (log(time_2014(i+1)) - log(time_2014(i)))^2;
end
vol(10) = sqrt(vol(10));

m = length(time_2015);
vol(11) = 0;
for i = 1:m-1
    vol(11) = vol(11) +(log(time_2015(i+1)) - log(time_2015(i)))^2;
end
vol(11) = sqrt(vol(11));

m = length(time_2016);
vol(12) = 0;
for i = 1:m-1
    vol(12) = vol(12)+(log(time_2016(i+1)) - log(time_2016(i)))^2;
end
vol(12) = sqrt(vol(12));

m = length(time_2017);
vol(13) = 0;
for i = 1:m-1
    vol(13) = vol(13) +(log(time_2017(i+1)) - log(time_2017(i)))^2;
end
vol(13) = sqrt(vol(13));

m = length(time_2018);
vol(14) = 0;
for i = 1:m-1
    vol(14) = vol(14) + (log(time_2018(i+1)) - log(time_2018(i)))^2;
end
vol(14) = sqrt(vol(14));


fprintf('2005, Yearly Return: %f, Yearly Volatility: %f \n',ret(1),vol(1));
fprintf('2006, Yearly Return: %f, Yearly Volatility: %f \n',ret(2),vol(2));
fprintf('2007, Yearly Return: %f, Yearly Volatility: %f \n',ret(3),vol(3));
fprintf('2008, Yearly Return: %f, Yearly Volatility: %f \n',ret(4),vol(4));
fprintf('2009, Yearly Return: %f, Yearly Volatility: %f \n',ret(5),vol(5));
fprintf('2010, Yearly Return: %f, Yearly Volatility: %f \n',ret(6),vol(6));
fprintf('2011, Yearly Return: %f, Yearly Volatility: %f \n',ret(7),vol(7));
fprintf('2012, Yearly Return: %f, Yearly Volatility: %f \n',ret(8),vol(8));
fprintf('2013, Yearly Return: %f, Yearly Volatility: %f \n',ret(9),vol(9));
fprintf('2014, Yearly Return: %f, Yearly Volatility: %f \n',ret(10),vol(10));
fprintf('2015, Yearly Return: %f, Yearly Volatility: %f \n',ret(11),vol(11));
fprintf('2016, Yearly Return: %f, Yearly Volatility: %f \n',ret(12),vol(12));
fprintf('2017, Yearly Return: %f, Yearly Volatility: %f \n',ret(13),vol(13));
fprintf('2018, Yearly Return: %f, Yearly Volatility: %f \n',ret(14),vol(14));


    

function [simup] = MCSimulation(ret,vol,time,close)

mu = ret(10);
theta = vol(10);
pos = time.Year == 2019;
n = length(pos(pos == 1));
close_2019 = close(pos);

simup = zeros(length(close_2019),5);
simup(1,:) = close_2019(1);

for iter = 1:5
    day = 1;
    for i = 0:0.5/118:1
        day = day+1;
        eps = normrnd(0,1);
        simup(day,iter) = exp(log(simup((day-1),iter))+(mu - theta^2/2)*(0.5/118)...
                      + theta * eps * sqrt(0.5/118));
    end
end






    
    

function resultplot(close_2019,simup)

m = length(close_2019);
plot(close_2019,'black');
hold on 
plot(simup(1:m,:));