adding hook for PATH solver

matlab/dynare_solve.m

@@ -184,6 +184,7 @@ elseif options.solve_algo == 3
         [x,info] = csolve(func,x,[],1e-6,500,varargin{:});
     end
 elseif options.solve_algo == 10
+    % LMMCP
     olmmcp = options.lmmcp;
 
     [x,fval,exitflag] = lmmcp(func,x,olmmcp.lb,olmmcp.ub,olmmcp,varargin{:});
@@ -192,7 +193,34 @@ elseif options.solve_algo == 10
     else
         info = 1;
     end
+elseif options.solve_algo == 11
+    % PATH linear mixed complementary problem
+    if ~exist('lcppath')
+        error(['PATH can''t be provided with Dynare. You need to install it ' ...
+               'yourself and add its location to Matlab/Octave path before ' ...
+               'running Dynare'])
+    end
+    if options.linear ~= 1
+        error(['solve_algo==11 (linear complementarity problem) requires ' ...
+               'linear  option in model instruction']);
+    end
+    olcppath = options.lcppath;
+    [junk,M] = func(x,varargin{:});
+    [x,mu] = pathlcp(fjac,olcppath.q,olcppath.lb,olcppath.ub,x,olcppath.A,olcppath.b,olcppath.t,olcppath.mu0);
+elseif options.solve_algo == 12
+    % PATH mixed complementary problem
+    % PATH linear mixed complementary problem
+    if ~exist('mcppath')
+        error(['PATH can''t be provided with Dynare. You need to install it ' ...
+               'yourself and add its location to Matlab/Octave path before ' ...
+               'running Dynare'])
+    end
+    omcppath = options.mcppath;
+    global mcp_data
+    mcp_data.func = func;
+    mcp_data.args = varargin;
+    [x,mu] = pathmcp(x,omcppath.lb,omcppath.ub,'mcp_func',omcppath.A,omcppath.b,omcppath.t,omcppath.mu0);
 else
-    error('DYNARE_SOLVE: option solve_algo must be one of [0,1,2,3,4,9,10]')
+    error('DYNARE_SOLVE: option solve_algo must be one of [0,1,2,3,4,9,10:12]')
 end
 

matlab/global_initialization.m

@@ -760,6 +760,18 @@ options_.convergence.geweke.geweke_interval=[0.2 0.5];
 % Options for lmmcp solver
 options_.lmmcp.status = 0;
 
+% Options for lcppath solver
+options_.lcppath.A = [];
+options_.lcppath.b = [];
+options_.lcppath.t = [];
+options_.lcppath.mu0 = [];
+
+% Options for mcppath solver
+options_.mcppath.A = [];
+options_.mcppath.b = [];
+options_.mcppath.t = [];
+options_.mcppath.mu0 = [];
+
 % initialize persistent variables in priordens()
 priordens([],[],[],[],[],[],1);
 % initialize persistent variables in dyn_first_order_solver()

matlab/mcp_func.m

@@ -0,0 +1,11 @@
+function [res,fjac,domer] = mcp_func(x,jacflag)
+global mcp_data
+
+domer = 0;
+if jacflag
+    [res,fjac] = mcp_data.func(x,mcp_data.args{:});
+else
+    res = mcp_data.func(x,mcp_data.args{:});
+    fjac = [];
+end
+disp(norm(res))

matlab/mcpath_function.m

@@ -0,0 +1,8 @@
+function [res,jac,domerr] = mcpath_function(func,z,jacflag,varargin)
+domerr = 0;
+
+if jacflag
+    [res,jac] = func(z,varargin{:});
+else
+    res = func(z,varargin{:});
+end

matlab/ms-sbvar/auto/sbvar_identification.el

@@ -0,0 +1,8 @@
+(TeX-add-style-hook "sbvar_identification"
+ (lambda ()
+    (TeX-run-style-hooks
+     "latex2e"
+     "art10"
+     "article"
+     "")))
+

matlab/perfect-foresight-models/perfect_foresight_solver_core.m

@@ -133,3 +133,5 @@ if nargout>1
         maxerror = max(max(abs(residuals)));
     end
 end
+
+

tests/path/lcp/mod1.mod

@@ -0,0 +1,63 @@
+var y pi rs rl rs_shadow rs_1;
+varexo e_y e_pi;
+
+parameters alpha_1 alpha_2 alpha_3 beta_1 beta_2 gamma_r gamma_pi gamma_y pi_bar rr_bar;
+
+
+alpha_1 = 0.7;
+alpha_2 = 0.99;
+alpha_3 = 0.1;
+beta_1 = 0.7;
+beta_2 = 0.08;
+gamma_r = 0.3;
+gamma_pi = 2;
+gamma_y = 1;
+pi_bar = 2.0;
+rr_bar = 1.0;
+
+model(linear,use_dll);
+    pi = pi_bar*(1-alpha_1-(1-alpha_1)*alpha_2) + alpha_1*pi(-1) + (1-alpha_1)*alpha_2*pi(+1) + alpha_3*y + e_pi;
+    y = beta_1*y(-1) + (1-beta_1)*y(+1) - beta_2*(rl-rr_bar - pi(+1)) + e_y;
+    rs_shadow = gamma_r*rs(-1) + (1-gamma_r)*(rr_bar + pi_bar + gamma_pi*(pi-pi_bar) + gamma_y*y);
+    [mcp = 'rs_shadow > 0']
+    rs_1 = rs_shadow;
+//    [mcp = 'pi > 1.4']
+    rs = rs_1;
+    rl = 400*(((1+rs/400)
+    @#for i in [1:39]
+        *(1+rs(+@{i})/400)
+    @#endfor
+    )^(1/40) - 1);
+end;
+
+steady_state_model;
+y = 0;
+pi = pi_bar;
+rs = rr_bar+pi_bar;
+rs_shadow = rs;
+rs_1 = rs;
+rl = rs;
+end;
+
+steady;
+
+check;
+
+shocks;
+var e_y;
+periods 1;
+//values -1.226;
+values -1;
+end;
+
+options_.solve_algo = 12;
+options_.mcp = 1;
+//options_.linear_approximation = 1;
+perfect_foresight_setup(periods=200);    
+
+perfect_foresight_solver(stack_solve_algo=7);
+
+rplot rs;
+//rplot rl;
+//rplot y;
+//rplot pi;

tests/path/lcp/mod1a.mod

@@ -0,0 +1,63 @@
+var y pi rs rl rs_shadow rs_1;
+varexo e_y e_pi;
+
+parameters alpha_1 alpha_2 alpha_3 beta_1 beta_2 gamma_r gamma_pi gamma_y pi_bar rr_bar;
+
+
+alpha_1 = 0.7;
+alpha_2 = 0.99;
+alpha_3 = 0.1;
+beta_1 = 0.7;
+beta_2 = 0.08;
+gamma_r = 0.3;
+gamma_pi = 2;
+gamma_y = 1;
+pi_bar = 2.0;
+rr_bar = 1.0;
+
+model(linear,use_dll);
+    pi = pi_bar*(1-alpha_1-(1-alpha_1)*alpha_2) + alpha_1*pi(-1) + (1-alpha_1)*alpha_2*pi(+1) + alpha_3*y + e_pi;
+    y = beta_1*y(-1) + (1-beta_1)*y(+1) - beta_2*(rl-rr_bar - pi(+1)) + e_y;
+    rs_shadow = gamma_r*rs(-1) + (1-gamma_r)*(rr_bar + pi_bar + gamma_pi*(pi-pi_bar) + gamma_y*y);
+//    [mcp = 'rs_shadow > 0']
+    rs_1 = rs_shadow;
+//    [mcp = 'pi > 1.4']
+    rs = rs_1;
+    rl = 400*(((1+rs/400)
+    @#for i in [1:39]
+        *(1+rs(+@{i})/400)
+    @#endfor
+    )^(1/40) - 1);
+end;
+
+steady_state_model;
+y = 0;
+pi = pi_bar;
+rs = rr_bar+pi_bar;
+rs_shadow = rs;
+rs_1 = rs;
+rl = rs;
+end;
+
+steady;
+
+check;
+
+shocks;
+var e_y;
+periods 1;
+values -1;
+end;
+
+//options_.solve_algo = 10;
+options_.mcp = 1;
+options_.linear_approximation = 1;
+perfect_foresight_setup(periods=200);    
+
+perfect_foresight_solver(stack_solve_algo=7);
+//perfect_foresight_solver;
+
+rplot rs;
+rplot rl;
+rplot y;
+rplot pi;