Script File to estimate log-periodic fit to USDZAR rand crash using MCMC
Data is USD/ZAR over 15000 days starting in 1960
Algorithm implemented in MCMC.M :
1. Choose initial point in parameter space 2. Generate parameter jump from jump-pdf 3. Accept/reject jump - always accept better likelihoods reject with uniform probability 4. Check burn-in and convergence/mixing
Function is:
From Sornette and Johansen, Quantitative Finance, 1, 452, 2001
Fitting Parameters are: (A, B, C, beta, t_c, w, phi)
dates are given in days since AD began but converted to approximate years using 365.25 d/yr
Bruce Bassett, Tim Gebbie
Contents
Clear the workspace
clear all;
clc;
Set initial conditions for parameter chain vector p
(A,B,C,BETA,TC,OMEGA,PHI)
para_str ={'A','B','C','\beta','t_c','\omega','\phi'}; epoch = '31-Jan-2000::31-Oct-2001'; p0 = [3.5,-3,2.27,0.35,2003,7,-14]; % initial parameters pu = [1,1,1,1,0,10,1]; % random weightings
Load the data
load data/workspace_zar.mat; % load file
Prepare the data
plot(ZAR); data = fts2mat(ZAR(epoch).USDZAR,1); ytilde = log(data(:,2)); % column data of ln(ZAR/USD) exchange rate (including NAN missing data) t = data(:,1)./ 365.25; % dates for the data (days since AD began) convert days to years chain_length = 3000;

Estimate the parameters
[chains,chi2]=mcmc(@logp,t,ytilde,p0,pu,chain_length);
Find the best theory
define theoretical prediction
best_theory = logp(t,chains(end,:));
Plot the best theory
figure;
plot(t, ytilde, t, best_theory);
% Check mixing/convergence using R-statistic (Verde et al)

Plot the chains
figure; for i=1:length(chains), subplot(4,2,i); plot(chains(:,i)); title(para_str{i}); end; figure for i=1:length(chains), subplot(4,2,i); hist(chains(:,i),40); title(para_str{i}); end;
Error using ==> evalin Index exceeds matrix dimensions.
