diff --git a/mex/sources/ms-sbvar/mex_ms_convert_free_parameters.cc b/mex/sources/ms-sbvar/mex_ms_convert_free_parameters.cc
index 2ea109462..e71943cd6 100644
--- a/mex/sources/ms-sbvar/mex_ms_convert_free_parameters.cc
+++ b/mex/sources/ms-sbvar/mex_ms_convert_free_parameters.cc
@@ -53,8 +53,8 @@ mexFunction(int nlhs, mxArray *plhs[],
   if (!mxIsDouble(prhs[1]))
     DYN_MEX_FUNC_ERR_MSG_TXT("Second argument is a vector of free parameters");
 
-  if (nlhs < 3)
-    DYN_MEX_FUNC_ERR_MSG_TXT("You must specify at least three output arguments [A0,Aplus,Zeta]");
+  if (nlhs < 4)
+    DYN_MEX_FUNC_ERR_MSG_TXT("You must specify at least four output arguments [err,A0,Aplus,Zeta]");
 
   // copy the string data from prhs[0] into a C string input_ buf.    */
   input_buf = mxArrayToString(prhs[0]);
@@ -86,17 +86,17 @@ mexFunction(int nlhs, mxArray *plhs[],
 
   // Ao (nstates x nvars x nvars)
   mwSize dims[3] = {nstates, nvars, nvars};
-  plhs[0] = mxCreateNumericArray(3, dims, mxDOUBLE_CLASS, mxREAL);
-  a0 = mxGetPr(plhs[0]);
+  plhs[1] = mxCreateNumericArray(3, dims, mxDOUBLE_CLASS, mxREAL);
+  a0 = mxGetPr(plhs[1]);
 
   // Zeta (nstates x nvars x nvars)
-  plhs[2] = mxCreateNumericArray(3, dims, mxDOUBLE_CLASS, mxREAL);
-  zeta = mxGetPr(plhs[2]);
+  plhs[3] = mxCreateNumericArray(3, dims, mxDOUBLE_CLASS, mxREAL);
+  zeta = mxGetPr(plhs[3]);
 
   // Aplus (nstates x (nlags*nvars +nconstant) x nvars)
   dims[1] = npre;
-  plhs[1] = mxCreateNumericArray(3, dims, mxDOUBLE_CLASS, mxREAL);
-  aplus = mxGetPr(plhs[1]);
+  plhs[2] = mxCreateNumericArray(3, dims, mxDOUBLE_CLASS, mxREAL);
+  aplus = mxGetPr(plhs[2]);
 
   // Grand Transition Matrix
   Q = mxCreateDoubleMatrix(nstates, nstates, mxREAL);
@@ -112,11 +112,11 @@ mexFunction(int nlhs, mxArray *plhs[],
     }
 
   // if they have passed 4 arguments then pass back Q
-  if (nlhs > 3)
-    plhs[3] = Q;
+  if (nlhs > 4)
+    plhs[4] = Q;
 
   mxFree(model);
-
+  plhs[0] = mxCreateDoubleScalar(0);
 }
 
 #endif
diff --git a/mex/sources/ms-sbvar/mex_ms_forecast.cc b/mex/sources/ms-sbvar/mex_ms_forecast.cc
index 9a7b3d5e6..678be6b95 100644
--- a/mex/sources/ms-sbvar/mex_ms_forecast.cc
+++ b/mex/sources/ms-sbvar/mex_ms_forecast.cc
@@ -54,8 +54,8 @@ mexFunction(int nlhs, mxArray *plhs[],
   // Check the left hand and right hand side arguments to make sure they conform
   if (mxIsChar(prhs[0]) != 1)
     DYN_MEX_FUNC_ERR_MSG_TXT("First argument has to be a string to the init_filename.");
-  if (nlhs != 1)
-    DYN_MEX_FUNC_ERR_MSG_TXT("You must specify one output argument");
+  if (nlhs != 2)
+    DYN_MEX_FUNC_ERR_MSG_TXT("You must specify two output arguments.");
 
   // copy the string data from prhs[0] into a C string input_ buf.    */
   input_buf = mxArrayToString(prhs[0]);
@@ -81,30 +81,30 @@ mexFunction(int nlhs, mxArray *plhs[],
       {
         // regimes x percentile x horizon x nvar
         mwSize dims[4] = {nstates, options->num_percentiles, options->horizon, nvars};
-        plhs[0] = mxCreateNumericArray(4, dims, mxDOUBLE_CLASS, mxREAL);
-        out_buf = mxGetPr(plhs[0]);
+        plhs[1] = mxCreateNumericArray(4, dims, mxDOUBLE_CLASS, mxREAL);
+        out_buf = mxGetPr(plhs[1]);
       }
     else
       {
         // regimes x horizon x (nvar*nvar)
         mwSize dims[3] = {nstates, options->horizon, nvars};
-        plhs[0] = mxCreateNumericArray(3, dims, mxDOUBLE_CLASS, mxREAL);
-        out_buf = mxGetPr(plhs[0]);
+        plhs[1] = mxCreateNumericArray(3, dims, mxDOUBLE_CLASS, mxREAL);
+        out_buf = mxGetPr(plhs[1]);
       }
   else
   if (options->num_percentiles > 1)
     {
       // percentile x horizon x (nvar*nvar)
       mwSize dims[3] = {options->num_percentiles, options->horizon, nvars};
-      plhs[0] = mxCreateNumericArray(3, dims, mxDOUBLE_CLASS, mxREAL);
-      out_buf = mxGetPr(plhs[0]);
+      plhs[1] = mxCreateNumericArray(3, dims, mxDOUBLE_CLASS, mxREAL);
+      out_buf = mxGetPr(plhs[1]);
     }
   else
     {
       // horizon x (nvar*nvar)
       mwSize dims[2] = {options->horizon, nvars};
-      plhs[0] = mxCreateNumericArray(2, dims, mxDOUBLE_CLASS, mxREAL);
-      out_buf = mxGetPr(plhs[0]);
+      plhs[1] = mxCreateNumericArray(2, dims, mxDOUBLE_CLASS, mxREAL);
+      out_buf = mxGetPr(plhs[1]);
     }
 
   // hard coding dates for now to be off
@@ -188,7 +188,7 @@ mexFunction(int nlhs, mxArray *plhs[],
 
   mxFree(p);
   mxFree(model);
-
+  plhs[0] = mxCreateDoubleScalar(0);
 }
 
 #endif
diff --git a/mex/sources/ms-sbvar/mex_ms_irf.cc b/mex/sources/ms-sbvar/mex_ms_irf.cc
index 3a2451264..ae607dd57 100644
--- a/mex/sources/ms-sbvar/mex_ms_irf.cc
+++ b/mex/sources/ms-sbvar/mex_ms_irf.cc
@@ -52,8 +52,8 @@ mexFunction(int nlhs, mxArray *plhs[],
   // Check the left hand and right hand side arguments to make sure they conform
   if (mxIsChar(prhs[0]) != 1)
     DYN_MEX_FUNC_ERR_MSG_TXT("First argument has to be a string to the init_filename.");
-  if (nlhs != 1)
-    DYN_MEX_FUNC_ERR_MSG_TXT("You must specify one output argument");
+  if (nlhs != 2)
+    DYN_MEX_FUNC_ERR_MSG_TXT("You must specify two output arguments.");
 
   // copy the string data from prhs[0] into a C string input_ buf.    */
   input_buf = mxArrayToString(prhs[0]);
@@ -78,30 +78,30 @@ mexFunction(int nlhs, mxArray *plhs[],
       {
         // regimes x percentile x horizon x (nvar*nvar)
         mwSize dims[4] = {nstates, options->num_percentiles, options->horizon, (nvars*nvars)};
-        plhs[0] = mxCreateNumericArray(4, dims, mxDOUBLE_CLASS, mxREAL);
-        out_buf = mxGetPr(plhs[0]);
+        plhs[1] = mxCreateNumericArray(4, dims, mxDOUBLE_CLASS, mxREAL);
+        out_buf = mxGetPr(plhs[1]);
       }
     else
       {
         // regimes x horizon x (nvar*nvar)
         mwSize dims[3] = {nstates, options->horizon, (nvars*nvars)};
-        plhs[0] = mxCreateNumericArray(3, dims, mxDOUBLE_CLASS, mxREAL);
-        out_buf = mxGetPr(plhs[0]);
+        plhs[1] = mxCreateNumericArray(3, dims, mxDOUBLE_CLASS, mxREAL);
+        out_buf = mxGetPr(plhs[1]);
       }
   else
   if (options->num_percentiles > 1)
     {
       // percentile x horizon x (nvar*nvar)
       mwSize dims[3] = {options->num_percentiles, options->horizon, (nvars*nvars)};
-      plhs[0] = mxCreateNumericArray(3, dims, mxDOUBLE_CLASS, mxREAL);
-      out_buf = mxGetPr(plhs[0]);
+      plhs[1] = mxCreateNumericArray(3, dims, mxDOUBLE_CLASS, mxREAL);
+      out_buf = mxGetPr(plhs[1]);
     }
   else
     {
       // horizon x (nvar*nvar)
       mwSize dims[2] = {options->horizon, (nvars*nvars)};
-      plhs[0] = mxCreateNumericArray(2, dims, mxDOUBLE_CLASS, mxREAL);
-      out_buf = mxGetPr(plhs[0]);
+      plhs[1] = mxCreateNumericArray(2, dims, mxDOUBLE_CLASS, mxREAL);
+      out_buf = mxGetPr(plhs[1]);
     }
 
   // Use filter probabilities?
@@ -188,7 +188,7 @@ mexFunction(int nlhs, mxArray *plhs[],
 
   mxFree(model);
   mxFree(p);
-
+  plhs[0] = mxCreateDoubleScalar(0);
 }
 
 #endif
diff --git a/mex/sources/ms-sbvar/mex_ms_variance_decomposition.cc b/mex/sources/ms-sbvar/mex_ms_variance_decomposition.cc
index 0a219859b..fdf0df038 100644
--- a/mex/sources/ms-sbvar/mex_ms_variance_decomposition.cc
+++ b/mex/sources/ms-sbvar/mex_ms_variance_decomposition.cc
@@ -51,8 +51,8 @@ mexFunction(int nlhs, mxArray *plhs[],
   /* Check the left hand and right hand side arguments to make sure they conform */
   if (mxIsChar(prhs[0]) != 1)
     DYN_MEX_FUNC_ERR_MSG_TXT("First argument has to be a string to the init_filename.");
-  if (nlhs != 1)
-    DYN_MEX_FUNC_ERR_MSG_TXT("You must specify one output argument");
+  if (nlhs != 2)
+    DYN_MEX_FUNC_ERR_MSG_TXT("You must specify two output arguments.");
 
   // copy the string data from prhs[0] into a C string input_ buf.    */
   input_buf = mxArrayToString(prhs[0]);
@@ -78,30 +78,30 @@ mexFunction(int nlhs, mxArray *plhs[],
       {
         /* regimes x percentile x horizon x (nvar*nvar) */
         mwSize dims[4] = {nstates, options->num_percentiles, options->horizon, (nvars*nvars)};
-        plhs[0] = mxCreateNumericArray(4, dims, mxDOUBLE_CLASS, mxREAL);
-        out_buf = mxGetPr(plhs[0]);
+        plhs[1] = mxCreateNumericArray(4, dims, mxDOUBLE_CLASS, mxREAL);
+        out_buf = mxGetPr(plhs[1]);
       }
     else
       {
         /* regimes x horizon x (nvar*nvar) */
         mwSize dims[3] = {nstates, options->horizon, (nvars*nvars)};
-        plhs[0] = mxCreateNumericArray(3, dims, mxDOUBLE_CLASS, mxREAL);
-        out_buf = mxGetPr(plhs[0]);
+        plhs[1] = mxCreateNumericArray(3, dims, mxDOUBLE_CLASS, mxREAL);
+        out_buf = mxGetPr(plhs[1]);
       }
   else
   if (options->num_percentiles > 1)
     {
       /* percentile x horizon x (nvar*nvar) */
       mwSize dims[3] = {options->num_percentiles, options->horizon, (nvars*nvars)};
-      plhs[0] = mxCreateNumericArray(3, dims, mxDOUBLE_CLASS, mxREAL);
-      out_buf = mxGetPr(plhs[0]);
+      plhs[1] = mxCreateNumericArray(3, dims, mxDOUBLE_CLASS, mxREAL);
+      out_buf = mxGetPr(plhs[1]);
     }
   else
     {
       /* horizon x (nvar*nvar) */
       mwSize dims[2] = {options->horizon, (nvars*nvars)};
-      plhs[0] = mxCreateNumericArray(2, dims, mxDOUBLE_CLASS, mxREAL);
-      out_buf = mxGetPr(plhs[0]);
+      plhs[1] = mxCreateNumericArray(2, dims, mxDOUBLE_CLASS, mxREAL);
+      out_buf = mxGetPr(plhs[1]);
     }
 
   /* Use filter probabilities? */
@@ -242,7 +242,7 @@ mexFunction(int nlhs, mxArray *plhs[],
 
   mxFree(p);
   mxFree(model);
-
+  plhs[0] = mxCreateDoubleScalar(0);
 }
 
 #endif