function check_contraction(delta,alpha);
% function check_contraction(delta,alpha);
%
% delta: rate of depeciation
% alpha: prameter in the production function
%
% function calculates both expressions in the contraction condition (8)
% and checks whether the inequality (8) holds
% Last update: April 30, 2001
%
% Author: Klaus Reiner Schenk-Hoppé 
% Email: klaus@iew.unizh.ch
% Web: www.iew.unizh.ch/home/klaus


format long

N = 10000;  % number of iterations
            % expected values are calculated according to the ergodic theorem
s = 0.25;

A = 0.95;
B = 0.95;
phi2 = 0.005;
rho2 = 0.0005;
rho1 = - rho2;
phi1 = - phi2;
nmin = 0.0075 + phi1/(1-A);
nmax = 0.0075 + phi2/(1-A);
smin = s * ( 1 + rho1/(1-B));
smax = s * ( 1 + rho2/(1-B));

eta(1) = 0.0;
xi(1) = 0.0;


%% generate the sample paths of the stochastic variables

rst = rand('state');

for t = 1:N,
  xi(t+1) = A * xi(t) + rho1 + (rho2-rho1)*rand;
  savings(t+1) = s*(1.0+xi(t+1));
  eta(t+1) = B * eta(t) + phi1 + (phi2-phi1)*rand;
  n(t) = 0.0075 + eta(t+1);
end

lhsresult = sum(log(1-delta + alpha./smin .* (delta+nmax) .* savings ))./length(savings);
rhsresult = sum(log(1+n))./length(n);

fprintf('Exp log (1-delta + alpha * s(omega)/s_min * (delta+n_max) ) = %f \n',lhsresult);
fprintf('Exp log (1+n(omega)) = %f \n',rhsresult);

if (lhsresult < rhsresult)
 fprintf('Contraction condition (8) is satisfied \n');
else 
 fprintf('Contraction condition (8) is NOT satisfied \n');
end