diff --git a/doc/dynare.texi b/doc/dynare.texi index ef146f4cf..8e24ac1fc 100644 --- a/doc/dynare.texi +++ b/doc/dynare.texi @@ -2984,8 +2984,14 @@ option, @pxref{Model declaration}) Trust-region algorithm on the entire model. @item 10 -Levenberg-Marquardt mixed complementarity problem (LMMCP) solver +Levenberg-Marquardt mixed compleproblem (LMMCP) solver (@cite{Kanzow and Petra 2004}) + +@item 11 +PATH 3.0 solver of @cite{Ferris and Munson (1999)}. Dynare only provides the interface +for using the solver. Due to licence restrictions, you have to download the solver yourself +from @url{http://pages.cs.wisc.edu/~ferris/path.html} and place it in Matlab's search path. + @end table @noindent @@ -13950,6 +13956,10 @@ Dynamic Equilibrium Economies: Linear versus Nonlinear Likelihood,'' Fernández-Villaverde, Jesús (2010): ``The econometrics of DSGE models,'' @i{SERIEs}, 1, 3--49 +@item +Ferris, Michael C. and Todd S. Munson (1999): ``Interfaces to PATH 3.0: Design, Implementation and Usage'', +@i{Computational Optimization and Applications}, 12(1), 207--227 + @item Geweke, John (1992): ``Evaluating the accuracy of sampling-based approaches to the calculation of posterior moments,'' in J.O. Berger, J.M. Bernardo, diff --git a/matlab/dynare_solve.m b/matlab/dynare_solve.m index f60b99d67..d8c1ee469 100644 --- a/matlab/dynare_solve.m +++ b/matlab/dynare_solve.m @@ -207,9 +207,12 @@ elseif options.solve_algo == 11 global mcp_data mcp_data.func = func; mcp_data.args = varargin; - - [x,fval,jac,mu,status] = pathmcp(x,omcppath.lb,omcppath.ub,'mcp_func',omcppath.A,omcppath.b,omcppath.t,omcppath.mu0); - info = ~status; + info=0; + try + [x,fval,jac,mu] = pathmcp(x,omcppath.lb,omcppath.ub,'mcp_func',omcppath.A,omcppath.b,omcppath.t,omcppath.mu0); + catch + info = 1; + end else error('DYNARE_SOLVE: option solve_algo must be one of [0,1,2,3,4,9,10,11]') end diff --git a/matlab/pathlcp.m b/matlab/pathlcp.m deleted file mode 100644 index 30ddc367e..000000000 --- a/matlab/pathlcp.m +++ /dev/null @@ -1,172 +0,0 @@ -function [z,mu,status] = pathlcp(M,q,l,u,z,A,b,t,mu) -% pathlcp(M,q,l,u,z,A,b,t,mu) -% -% Solve the standard linear complementarity problem using PATH: -% z >= 0, Mz + q >= 0, z'*(Mz + q) = 0 -% -% Required input: -% M(n,n) - matrix -% q(n) - vector -% -% Output: -% z(n) - solution -% mu(m) - multipliers (if polyhedral constraints are present) -% -% Optional input: -% l(n) - lower bounds default: zero -% u(n) - upper bounds default: infinity -% z(n) - starting point default: zero -% A(m,n) - polyhedral constraint matrix default: empty -% b(m) - polyhedral right-hand side default: empty -% t(m) - type of polyhedral constraint default: 1 -% < 0: less than or equal -% 0: equation -% > 0: greater than or equal -% mu(m) - starting value for multipliers default: zero -% -% The optional lower and upper bounds are used to define a linear mixed -% complementarity problem (box constrained variational inequality). -% l <= z <= u -% where l_i < z_i < u_i => (Mz + q)_i = 0 -% l_i = z => (Mz + q)_i >= 0 -% u_i = z => (Mz + q)_i <= 0 -% -% The optional constraints are used to define a polyhedrally constrained -% variational inequality. These are transformed internally to a standard -% mixed complementarity problem. The polyhedral constraints are of the -% form -% Ax ? b -% where ? can be <=, =, or >= depending on the type specified for each -% constraint. - -Big = 1e20; - -if (nargin < 2) - error('two input arguments required for lcp(M, q)'); -end - -if (~issparse(M)) - M = sparse(M); % Make sure M is sparse -end -q = full(q(:)); % Make sure q is a column vector - -[mm,mn] = size(M); % Get the size of the inputs -n = length(q); - -if (mm ~= mn | mm ~= n) - error('dimensions of M and q must match'); -end - -if (n == 0) - error('empty model'); -end - -if (nargin < 3 | isempty(l)) - l = zeros(n,1); -end - -if (nargin < 4 | isempty(u)) - u = Big*ones(n,1); -end - -if (nargin < 5 | isempty(z)) - z = zeros(n,1); -end - -z = full(z(:)); l = full(l(:)); u = full(u(:)); -if (length(z) ~= n | length(l) ~= n | length(u) ~= n) - error('Input arguments are of incompatible sizes'); -end - -l = max(l,-Big*ones(n,1)); -u = min(u,Big*ones(n,1)); -z = min(max(z,l),u); - -m = 0; -if (nargin > 5) - if (nargin < 7) - error('Polyhedral constraints require A and b'); - end - - if (~issparse(A)) - A = sparse(A); - end - b = full(b(:)); - - m = length(b); - - if (m > 0) - - [am, an] = size(A); - - if (am ~= m | an ~= n) - error('Polyhedral constraints of incompatible sizes'); - end - - if (nargin < 8 | isempty(t)) - t = ones(m,1); - end - - if (nargin < 9 | isempty(mu)) - mu = zeros(m,1); - end - - t = full(t(:)); mu = full(mu(:)); - if (length(t) ~= m | length(mu) ~= m) - error('Polyhedral input arguments are of incompatible sizes'); - end - - l_p = -Big*ones(m,1); - u_p = Big*ones(m,1); - - idx = find(t > 0); - if (length(idx) > 0) - l_p(idx) = zeros(length(idx),1); - end - - idx = find(t < 0); - if (length(idx) > 0) - u_p(idx) = zeros(length(idx),1); - end - - mu = min(max(mu,l_p),u_p); - - M = [M -A'; A sparse(m,m)]; - q = [q; -b]; - - z = [z; mu]; - l = [l; l_p]; - u = [u; u_p]; - else - if (nargin >= 9 & ~isempty(mu)) - error('No polyhedral constraints -- multipliers set.'); - end - - if (nargin >= 8 & ~isempty(t)) - error('No polyhedral constraints -- equation types set.'); - end - end -end - -idx = find(l > u); -if length(idx) > 0 - error('Bounds infeasible.'); -end - -nnzJ = nnz(M); - -[status, ttime] = lcppath(n+m, nnzJ, z, l, u, M, q); - -%if (status ~= 1) -% status, -% error('Path fails to solve problem'); -%end - -mu = []; -if (m > 0) - mu = z(n+1:n+m); - z = z(1:n); -end - -return; - diff --git a/matlab/pathmcp.m b/matlab/pathmcp.m deleted file mode 100644 index 1db4dccb1..000000000 --- a/matlab/pathmcp.m +++ /dev/null @@ -1,197 +0,0 @@ -function [z,f,J,mu,status] = pathmcp(z,l,u,cpfj,A,b,t,mu) -% pathmcp(z,l,u,cpfj,A,b,t,mu) -% -% Solve a polyhedrally constrained variational inequality using PATH -% -% Calling syntax: [z,f,J] = pathmcp(z,l,u,cpfunjac,A,b,t,mu) -% -% Input: -% z - starting point -% l - lower bounds on z -% u - upper bounds on z -% -% cpfunjac - the name of the m-file for evaluating the function F and its -% Jacobian J (without .m-extension). -% -% The following m-file must be supplied (where default name is -% 'mcp_funjac.m' unless stated otherwise in the variable cpfunjac). -% -% 'mcp_funjac.m' contains function [f,J,domerr]=cpfunjac(z,jacflag) -% that computes the function F and if jacflag=1 the sparse -% Jacobian J at the point z. domerr returns the number of domain -% violations. -% -% A - constraint matrix -% b - right hand side of the constraints -% t - types of the constraints -% <0 : less than or equal -% =0 : equal to -% >0 : greater than or equal -% -% We have Ax ? b, ? is the type of constraint -% -% Output: -% z - solution -% mu - multipliers on the constraints -% f - function evaluation at the solution -% J - jacobian evaluation at the solution - - -Big = 1e20; - -if (nargin < 1) - error('one input arguments required for mcp(z)'); -end - -z = full(z(:)); -n = length(z); - -if (n == 0) - error('empty model'); -end - -if (nargin < 2 | isempty(l)) - l = zeros(n,1); -end - -if (nargin < 3 | isempty(u)) - u = Big*ones(n,1); -end - -l = full(l(:)); u = full(u(:)); -if (length(l) ~= n | length(u) ~= n) - error('Input arguments are of incompatible sizes'); -end - -l = max(l,-Big*ones(n,1)); -u = min(u,Big*ones(n,1)); -z = min(max(z,l),u); - -if (nargin < 4 | isempty(cpfj)) - cpfj = 'mcp_funjac'; -end - -m = 0; -mu = []; -l_p = []; -u_p = []; - -if (nargin > 4) - if (nargin < 6) - error('Polyhedral constraints require A and b'); - end - - if (~issparse(A)) - A = sparse(A); - end - b = full(b(:)); - - m = length(b); - - if (m > 0) - - [am, an] = size(A); - - if (am ~= m | an ~= n) - error('Polyhedral constraints of incompatible sizes'); - end - - if (nargin < 7 | isempty(t)) - t = ones(m,1); - end - - if (nargin < 8 | isempty(mu)) - mu = zeros(m,1); - end - - t = full(t(:)); mu = full(mu(:)); - if (length(t) ~= m | length(mu) ~= m) - error('Polyhedral input arguments are of incompatible sizes'); - end - - l_p = -Big*ones(m,1); - u_p = Big*ones(m,1); - - idx = find(t > 0); - if (length(idx) > 0) - l_p(idx) = zeros(length(idx),1); - end - - idx = find(t < 0); - if (length(idx) > 0) - u_p(idx) = zeros(length(idx),1); - end - - mu = min(max(mu,l_p),u_p); - else - if (nargin >= 8 & ~isempty(mu)) - error('No polyhedral constraints -- multipliers set.'); - end - - if (nargin >= 7 & ~isempty(t)) - error('No polyhedral constraints -- equation types set.'); - end - end -else - A = []; -end - -% this is a fix, nnz may be bigger than this -[f,J,domerr] = feval(cpfj,z+1e-5*ones(size(z))+1e-5*abs(z),1); - -if (domerr > 0) - [f,J,domerr] = feval(cpfj,z,1); -end - -if (domerr > 0) - error([cpfj ' not defined at starting point']); -end - -if ~issparse(J) - error([cpfj ' must return a sparse Jacobian']); -end - -nnzJ = nzmax(J); - -row = n + m; -ele = nnzJ + 2*nzmax(A); - -init = [z; mu]; -low = [l; l_p]; -upp = [u; u_p]; - -if m > 0 - global mcp_vifunc; - global mcp_viconn; - global mcp_viconm; - global mcp_viconA; - global mcp_viconb; - - mcp_vifunc = cpfj; - mcp_viconn = n; - mcp_viconm = m; - mcp_viconA = A; - mcp_viconb = b; - - [status, ttime, f, J] = mcppath(row, ele, init, low, upp, 'mcp_vifunjac'); -else - [status, ttime, f, J] = mcppath(row, ele, init, low, upp, cpfj); -end - -%if (status ~= 1) -% status, -% error('Path fails to solve problem'); -%end - -mu = []; -z = init; - -if m > 0 - mu = init(n+1:n+m); - z = init(1:n); - - J = J(1:n,1:n); - f = f(1:n) + A'*mu; -end - -return;