diff --git a/matlab/default_option_values.m b/matlab/default_option_values.m
index 8cc2428f291badca82c8ccc03131261d6bea1b30..db52fe5902e12c9ac1300714a65c2dc983f2026d 100644
--- a/matlab/default_option_values.m
+++ b/matlab/default_option_values.m
@@ -72,6 +72,7 @@ options_.huge_number = 1e7;
 options_.threads.kronecker.sparse_hessian_times_B_kronecker_C = num_procs;
 options_.threads.local_state_space_iteration_2 = 1;
 options_.threads.perfect_foresight_problem = num_procs;
+options_.threads.k_order_perturbation = max(1, num_procs/2);
 
 % steady state
 options_.jacobian_flag = true;
diff --git a/matlab/set_dynare_threads.m b/matlab/set_dynare_threads.m
index ca7827cb90c79829e7f419292e0b530107ee70b2..300b00c970e0c2187d47ee9139eee4ec0b5e099a 100644
--- a/matlab/set_dynare_threads.m
+++ b/matlab/set_dynare_threads.m
@@ -42,6 +42,8 @@ switch mexname
     options_.threads.local_state_space_iteration_2 = n;
   case 'perfect_foresight_problem'
     options_.threads.perfect_foresight_problem = n;
+  case 'k_order_perturbation'
+    options_.threads.k_order_perturbation = n;
   otherwise
     message = [ mexname ' is not a known parallel mex file.' ];
     message_id  = 'Dynare:Threads:UnknownParallelMex';
diff --git a/mex/sources/k_order_perturbation/k_order_perturbation.cc b/mex/sources/k_order_perturbation/k_order_perturbation.cc
index 32feff67ac4a11b9b1ee44adbfc96a1770585deb..03578a97a4abbd48c92642bd8d657ff7f61e8057 100644
--- a/mex/sources/k_order_perturbation/k_order_perturbation.cc
+++ b/mex/sources/k_order_perturbation/k_order_perturbation.cc
@@ -113,6 +113,14 @@ extern "C" {
     if (qz_criterium_mx && mxIsScalar(qz_criterium_mx) && mxIsNumeric(qz_criterium_mx))
       qz_criterium = mxGetScalar(qz_criterium_mx);
 
+    const mxArray *threads_mx = mxGetField(options_mx, 0, "threads");
+    if (!threads_mx)
+      DYN_MEX_FUNC_ERR_MSG_TXT("Can't find field options_.threads");
+    const mxArray *num_threads_mx = mxGetField(threads_mx, 0, "k_order_perturbation");
+    if (!(num_threads_mx && mxIsScalar(num_threads_mx) && mxIsNumeric(num_threads_mx)))
+      DYN_MEX_FUNC_ERR_MSG_TXT("options_.threads.k_order_perturbation be a numeric scalar");
+    int num_threads = static_cast<int>(mxGetScalar(num_threads_mx));
+
     // Extract various fields from M_
     const mxArray *fname_mx = mxGetField(M_mx, 0, "fname");
     if (!(fname_mx && mxIsChar(fname_mx) && mxGetM(fname_mx) == 1))
@@ -200,6 +208,9 @@ extern "C" {
         // intiate tensor library
         TLStatic::init(kOrder, nStat+2*nPred+3*nBoth+2*nForw+nExog);
 
+        // Set number of parallel threads
+        sthread::detach_thread_group::max_parallel_threads = num_threads;
+
         // make KordpDynare object
         KordpDynare dynare(endoNames, exoNames, nExog, nPar,
                            ySteady, vCov, modParams, nStat, nPred, nForw, nBoth,