% This code illustrates Proposition 6 in Alvarez, Lebihan and Lippi (2014)
% "Small and large price changes and the propagation of monetary shocks"
%
% It computes the cumulated real effect of monetary policy as given by
% Mprime = \partial M(0) / \partial delta , and compares it with
% Kurt(Dp)/ ( 6*N(Dp) ) .
% Computes Mprime by evaluating the analytic expression given in the paper.
%
% The code used several functions, each one with a name fct_*.m where *
% is the rest of the name. See list of functions below
clear all
close all
%%% Input parameters 4 parameters
lambda = 0.5; % arrival rate of free adjustment, a strictly positive number
n = 10 ; % number of products, an integer equal or larger than one
ybar = 1/10 ; % threshold of the norm square of the n price gaps
sigma2 = 0.1^2; % variance of shocks, a strictly positive number
% to check the equality eqn 30 in Lemma 5 only the ratio (lambda * ybar / sigma2 ) matters
% We normalize epsilon to one since (1/epsilon) multiple every final expression
epsilon= 1 ; % curvature of utility; CRRA parameter
Nterms=30; % number of terms in the sumes to be included in different power series
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% start of the code %%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% To check eq 30 of lemma 5 of Proposition 6 one can use one function fct_rhs_lhs.m
% [lhs rhs rhs1 rhs2] = fct_rhs_lhs(lambda,ybar,n,sigma2,Nterms);
% the LHS of Lemma 5 is sum(lhs)
% the RHS of Lemma 5 is sum(rhs)
disp(' ') ;
disp( ' Parameters: ')
disp( [ ' n = ' ,num2str(n,2), ', lambda = ' num2str(lambda,2), ', sigma = ', num2str(sigma2^.5,2), ...
', bar y = ',num2str(ybar,2) ] );
disp(' ') ;
disp(' LHS and RHS of equation 30 equivalent to lambda Kur/Na 6 = lambda int [T+Tp y2/n] f dy')
disp(' i.e. LHS and RHS of equation 30 in Lemma 5 in the proof of Proposition 6 :')
[lhs rhs rhs1 rhs2] = fct_rhs_lhs(lambda,ybar,n,sigma2,Nterms);
disp([' LHS of equation 30 = ', num2str(sum(lhs),4)]); % LHS of equation 25
disp([' RHS of equation 30 = ', num2str(sum(rhs),4)]); % RHS of equation 25
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% For completness we also compute several statistics in alternative ways
% compute N(Dp) using Na =1/T(0)
Na = 1/fct_Ty(0,lambda,n,ybar,sigma2,Nterms) ;
% compute kurtosis by using the reprsetation of f and the forumla for moments of Dp
% in term of weighte sum of Bessel fcts
nu = abs(n/2-1); % order of Bessel functions
% CI and CK are weight of bessel function I and bessel function K
ratio = besseli(nu,2*sqrt(ybar*lambda/(2*sigma2)))/besselk(nu,2*sqrt(ybar*lambda/(2*sigma2)));
CI = 1/( quadgk(@(y) func_besseli_int(y,sigma2^.5,lambda,n), 0, ybar) ...
- ratio* quadgk(@(y) func_besselk_int(y,sigma2^.5,lambda,n), 0, ybar) ) ;
CK = -ratio*CI;
% computes integral in 4rth moment
mom_y2 = quadgk(@(y) y.^2.*CI.*func_besseli_int(y,sigma2^.5,lambda,n) + ...
y.^2.*CK.*func_besselk_int(y,sigma2^.5,lambda,n), 0, ybar);
% computes integrals for second moment
mom_y1 = quadgk(@(y) y.*CI.*func_besseli_int(y,sigma2^.5,lambda,n) + ...
y.*CK.*func_besselk_int(y,sigma2^.5,lambda,n), 0, ybar);
% compute kurtosis
kurt = 3* n/(2+n) * ( (Na-lambda)/Na * ybar^2 + lambda/Na* mom_y2 ) ...
./ ( (Na-lambda)/Na * ybar + lambda/Na* mom_y1 ).^2 ;
% Now compute M'(0) using direct integration (as in Lemma 3)
Mp0 = quad( @(z) (fct_Ty(z,lambda,n+2,ybar,sigma2,Nterms)+ 2/n.*z.*fct_Tpy(z,lambda,n+2,ybar,sigma2,Nterms))...
.*fct_f_multi(z,sigma2^.5,lambda,n,CI,CK),0,ybar)/epsilon ;
% alternative computaiton of M'(0) using a power function representation of f
Mp0_alt = quad( @(z) (1/lambda)*fct_lam_TpTpry2n(z,lambda,ybar,n,sigma2^.5,Nterms) ...
.* fct_f_poly_odd_n(z,lambda,n,ybar,sigma2,Nterms), 0,ybar ) / epsilon ;
% computes Kurtosis over N using power function representation:
kurt_over_Na = fct_kurt_over_N(lambda,ybar,n,sigma2^.5,Nterms);
% displays output in screen
disp(' ') ;
disp(' Other statistics' )
disp(' ') ;
disp([' Fraction free adjustments \ell = ',num2str(lambda/Na,2 ) ])
disp([' kurtosis price changes = ', num2str(kurt,2)]) ;
disp([' number price changes = ', num2str(Na,2)]) ;
disp(' ') ;
disp([ ' Mp(0) (derivative of M(d) evaluated at d=0 is ',num2str(Mp0,4 ) ])
disp([ ' the ratio kurt/(epsilon*6*Na) is ',num2str(kurt/(epsilon*6*Na),4 ) ])
disp([ ' the ratio kurt/(6) is ',num2str(kurt/(6),4 ),' and Mp(0)* epsilon * Na is ',num2str(Mp0*Na*epsilon,4 ) ])
disp(' ') ;
disp([ ' using power series: Mp(0) ( derivative of M(d) evaluated at d=0 ) is ',num2str(Mp0_alt,4 ) ]) ;
disp([ ' using power series: the ratio kurt/(epsilon*6*Na) is ',num2str(kurt_over_Na/6,4 ) ]) ;
disp(' ');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% List of functions used in the code:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% besseli, besselk, quadgk (these are public domain)
% created by us:
% func_besseli_int, func_besselk_int, fct_Ty, fct_Tpy, fct_kurt_over_N, fct_rhs_lhs
% fct_f_poly_odd_n, fct_int_yj_f_poly_odd_n, fct_lam_TpTpry2n, fct_f_multi