%probe.m
clear all 
format compact

GP = 100; %number of grid points

PR = 1/100; %probability of retrieving probe to local peak 

get_d;
%The last 4 elements of the second column of D  were obtained as std(Y)'*0.25;  This implies that the upper bound of the prior distribution of the standard  deviation of measurement errors is 25 percent of the standard deviation of the observables.  Equivalently, measurement errors are allowed to explain at most 6.25 percent of the variance of the observable. 

LB = D(:,1); %lower bound
UB = D(:,2);%upper bound

param0 = (LB+UB)/2;
%load probe.mat mxpar
%param0 = mxpar;


step = (UB-LB)/(GP-1);

mxpar = param0;
MXPAR =param0; %chain of assoc. params

%Number of estimated parameters
np = length(param0);

%load data
Y1 = xlsread('mexico_data.xls');
%columns: date, growth rate of output, growth rate of consumption, growth rate of investment, trade-balance-to-gdp ratio.
%Note for the authors: mexico_data.xls was created by running data_level2growth.m in
%c:\data\uribe\book\noise_shocks\mfiles
T = size(Y1,1);  %number of observations
%Demean the data and remove date column
Y = Y1(:,2:end) - ones(T,1)*mean(Y1(:,2:end));

pick_variables = [1 2 3 4]; %picks the elements of the policy funciton (rows of gx) corresponding to the observables.

%Initial Prior likelihood
logprior0 = log(prod( duniform(param0,D(:,1),D(:,2))));
 
%Calculate the initial log likelihood of the model for param0
loglikelihood0 =  log_likelihood(param0,pick_variables,Y);
 
%calculate initial posterior density value
mxlp =  loglikelihood0 + logprior0;
MXLP = mxlp; %chain of max logposteriors

save probe.mat

param = mxpar;
nmx = 1; %number of maxima found

while 2>1

%Draw parameter value
 rand(np,1);
(ans>2/3)-(ans<1/3);
 param = param +ans .* rand(np,1) .* step;
param = min(param,UB);
param = max(param,LB);

if rand<PR
param = mxpar;
end


lastwarn('');
loglikelihood =  log_likelihood(param,pick_variables,Y);

[a,b] = lastwarn;
if ~strcmp(b,'')
loglikelihood = -inf
end

logposterior = loglikelihood + logprior0;

%Pick param if higher likelihoods values
if logposterior>mxlp
mxlp =  logposterior;
MXLP = [MXLP mxlp];
mxl = loglikelihood
mxpar = param(:);
MXPAR = [MXPAR mxpar];
nmx = nmx+1
PR =  1/  (100 + 900*(rand>0.9));
save probe.mat 
end %if logposterior>mxlp
end %while 2>1