Code covered by the BSD License

# Mining Economics with MATLAB

### David Willingham (view profile)

03 Dec 2012 (Updated )

Simulates a distribution of Net Present Values of a mine via a forecasting model for Iron Ore Prices

ironoreprice_script.m
%% Simulating Iron Ore Prices
% This example demonstrates calibrating an geomtric brownian motion &
% mean reverting stochastic model from historical data of iron ore prices.
% The model is then used to simulate the prices into the future using
% the Stochastic Differential Equation Simulation engine in Econometrics
% Toolbox.
%
%% Import Historical Data
% Data is saved as a excel file based on historical iron ore prices from
% May 10, 2012 to June 22, 2012
% http://www.indexmundi.com/commodities/?commodity=iron-ore&months=240&currency=aud

% Database Version
importirondb
dates = data.Month;
Date = datenum(dates,'yyyy-mm-dd HH:MM:SS.FFF');
S = data.Price;

% Excel
% dates = txt(3:end,1);
% Date = datenum(dates,'dd/mm/yyyy');
% S = data(:,1);
SR = S(2:end)./S(1:end-1) - 1;
%% The Models
% Return Model (gbm)
% The returns based model used for simulation Iron Orce Prices is geometric
% brownian motion.  It can be described as:
%
% $$dX_t = \mu(t)X_t dt + D(t,X_t)V(t)dW_t$$
%
% $$\mu = \textrm{ Instantaneous Rate of Return}$$
%
% $$\sigma = \textrm { Volatility Rate}$$
%
% Mean Reversion Model (hwv)
% The mean reversion model used for simulating the Iron Ore prices is an Ornstein-Uhlenbeck
% brownian motion with mean reverting drift. This model is fit to the prices
% iron ore prices. The discrete-time equation of this model can be written as,
%
% $$\Delta x_t = \alpha(\mu - x_t)\Delta t + \sigma dz_t, \textrm{ where } dz_t \sim N(0, \sqrt\Delta_t)$$
%
% $$\alpha = \textrm{ Mean Reversion Rate}$$
%
% $$\mu = \textrm{ Mean level}$$
%
% $$\sigma = \textrm { Volatility }$$

%% Calibrate Parameters

modelsel = 'gbm'; %hwv - mean reverting, %gbm - based on returns
startdate = '2010-08-01';
i = find(datenum(startdate,'yyyy-mm-dd') == Date);
S2 = S(i:end);
Date2 = Date(i:end);
SR2 = SR(i+1:end);

switch modelsel
case 'gbm'
dt = 1;
sigma = std(SR2);
mu = mean(SR2);
model = gbm(mu,sigma,'StartState', S(end));
case 'hwv'
% hwv - mean reversion model
% The reversion rate and mean level can be calculated from the coefficients
% of a linear fit between the prices and their first difference scaled
% by the time interval parameter. All quantities are specified on an annual
% scale.
regressors = [ones(length(S2) - 1, 1) S2(1:end-1)];
[coefficients, intervals, residuals] = regress(diff(S2), regressors);
dt    = 1;
speed = -coefficients(2)/dt;
level = -coefficients(1)/coefficients(2);
sigma =  std(residuals)/sqrt(dt);
model = hwv(speed,level,sigma,  'StartState', S(end));
end

%% Monte-Carlo Simulation
% The model defined above can be simulated with the simulate method of the
% SDE object to generate multiple log price paths. These are exponentiated
% to compute the simulated natural gas prices. The plot below shows 1000
% paths simulated daily, 8 years into the future.

NTrials = 1000;
years = 8;
ann = 12;
NSteps  = years*ann;

tic
Xsim = simulate(model, NSteps, 'NTrials', NTrials, 'DeltaTime', dt);
Xsim = squeeze(Xsim); % Remove redundant dimension
toc

Ssim = Xsim;
SYear = Ssim(13:12:end,:);

% Visualize first 20 prices of 1000 paths
plot(Date, S, Date(end)+(0:30:(years*ann+1)*30-1), Ssim(:,1:100));
datetick; xlabel('Date'); ylabel('Iron Ore Price $US/tonne'); title('Historical & Simulated Iron Ore Prices') %% Save Model % The calibrated model is saved in a MAT-file for later use. save ironoremodel model dt %% Visual Analysis of Simulated Price Paths % Instead of plotting a number of paths at once, we can plot longer single paths % against the observed historical data to visually validate the simulated % paths. This can serve as a final sanity check. path = 15; figure plot(Date, S, 'b', Date(end)+(0:30:NSteps*30), Ssim(:,path), 'r'); title(['Historical & Simulated Prices, Path ' int2str(path)]); datetick('x','keeplimits'); xlabel('Date'); ylabel('Iron Ore Price$US/tonne');
%% Histogram results of price paths at the end of each year

figure

for i = 1:years
subplot(years/2,2,i)
SsimM = mean(Ssim(i*12+1,:));
hist(Ssim(i*12,:),20)
hold on
hm = plot(SsimM,0,'ro');
title(['Year: ',num2str(i),' of Iron Ore Prices'])
xlabel(['Iron Ore Price: Mean =',num2str(round(100*(SsimM))/100)])
zlabel('Simulation Count')
end