diff --git a/matlab/ols/dyn_ols.m b/matlab/ols/dyn_ols.m
index 411df6c8f..2f5a9379a 100644
--- a/matlab/ols/dyn_ols.m
+++ b/matlab/ols/dyn_ols.m
@@ -51,6 +51,7 @@ end
 
 %% Get Equation(s)
 jsonmodel = loadjson(jsonfile);
+ast       = jsonmodel.abstract_syntax_tree;
 jsonmodel = jsonmodel.model;
 
 if nargin == 1
@@ -61,103 +62,137 @@ elseif nargin == 2
         (size(fitted_names_dict, 2) == 2 || size(fitted_names_dict, 2) == 3)), ...
         'dyn_ols: the second argument must be an Nx2 or Nx3 cell array');
 elseif nargin == 3
+    ast = getEquationsByTags(ast, 'name', eqtags);
     jsonmodel = getEquationsByTags(jsonmodel, 'name', eqtags);
 end
 
 %% Estimation
-M_endo_exo_names_trim = [M_.endo_names; M_.exo_names];
-[junk, idxs] = sort(cellfun(@length, M_endo_exo_names_trim), 'descend');
-regex = strjoin(M_endo_exo_names_trim(idxs), '|');
-mathops = '[\+\*\^\-\/\(\)]';
-st = dbstack(1);
-varargout = cell(1, 1);
-called_from_olsgibbs = false;
-if strcmp(st(1).name, 'olsgibbs')
-    varargout = cell(1, 4);
-    called_from_olsgibbs = true;
-end
-for i = 1:length(jsonmodel)
-    %% Construct regression matrices
-    rhs_ = strsplit(jsonmodel{i}.rhs, {'+','-','*','/','^','log(','exp(','(',')'});
-    rhs_(cellfun(@(x) all(isstrprop(x, 'digit')), rhs_)) = [];
-    vnames = setdiff(rhs_, M_.param_names);
-    if ~isempty(regexp(jsonmodel{i}.rhs, ...
-            ['(' strjoin(vnames, '\\(\\d+\\)|') '\\(\\d+\\))'], ...
-            'once'))
-        error(['dyn_ols: you cannot have leads in equation on line ' ...
-            jsonmodel{i}.line ': ' jsonmodel{i}.lhs ' = ' jsonmodel{i}.rhs]);
+
+for i = 1:length(ast)
+    if ~strcmp(ast{i}.AST.node_type, 'BinaryOpNode') ...
+            || ~strcmp(ast{i}.AST.op, '=')
+        ols_error('expecting equation with equal sign', line);
     end
 
-    pnames = intersect(rhs_, M_.param_names);
-    vnames = cell(1, length(pnames));
-    splitstrings = cell(length(pnames), 1);
+    % Check LHS
+    if ~strcmp(ast{i}.AST.arg1.node_type, 'VariableNode') ...
+            && ~strcmp(ast{i}.AST.arg1.node_type, 'UnaryOpNode')
+        ols_error('expecting Variable or UnaryOp on LHS', line);
+    else
+        if strcmp(ast{i}.AST.arg1.node_type, 'VariableNode') ...
+                && (ast{i}.AST.arg1.lag ~= 0 ...
+                || ~any(strcmp(ds.name, ast{i}.AST.arg1.name)))
+            ols_error('the lhs of the equation must be an Variable or UnaryOp with lag == 0 that exists in the dataset', line);
+        end
+        if strcmp(ast{i}.AST.arg1.node_type, 'UnaryOpNode') ...
+                && (ast{i}.AST.arg1.arg.lag ~= 0 ...
+                || ~any(strcmp(ds.name, ast{i}.AST.arg1.arg.name)))
+            ols_error('the lhs of the equation must be an Variable or UnaryOp with lag == 0 that exists in the dataset', line);
+        end
+    end
+
+    % Set LHS (Y)
+    lhssub = dseries();
+    Y = evalNode(ds, ast{i}.AST.arg1, jsonmodel{i}.line, dseries());
+
+    % Set RHS (X)
+    plus_node = ast{i}.AST.arg2;
+    last_node_to_parse = [];
+    residual = '';
     X = dseries();
-    for j = 1:length(pnames)
-        createdvar = false;
-        pregex = [...
-            mathops pnames{j} mathops ...
-            '|^' pnames{j} mathops ...
-            '|' mathops pnames{j} '$' ...
-            ];
-        [startidx, endidx] = regexp(jsonmodel{i}.rhs, pregex, 'start', 'end');
-        assert(length(startidx) == 1);
-        if jsonmodel{i}.rhs(startidx) == '*' && jsonmodel{i}.rhs(endidx) == '*'
-            vnamesl = getStrMoveLeft(jsonmodel{i}.rhs(1:startidx-1));
-            vnamesr = getStrMoveRight(jsonmodel{i}.rhs(endidx+1:end));
-            vnames{j} = [vnamesl '*' vnamesr];
-            splitstrings{j} = [vnamesl '*' pnames{j} '*' vnamesr];
-        elseif jsonmodel{i}.rhs(startidx) == '*'
-            vnames{j} = getStrMoveLeft(jsonmodel{i}.rhs(1:startidx-1));
-            splitstrings{j} = [vnames{j} '*' pnames{j}];
-        elseif jsonmodel{i}.rhs(endidx) == '*'
-            vnames{j} = getStrMoveRight(jsonmodel{i}.rhs(endidx+1:end));
-            splitstrings{j} = [pnames{j} '*' vnames{j}];
-            if jsonmodel{i}.rhs(startidx) == '-'
-                vnames{j} = ['-' vnames{j}];
-                splitstrings{j} = ['-' splitstrings{j}];
-            end
-        elseif jsonmodel{i}.rhs(startidx) == '+' ...
-                || jsonmodel{i}.rhs(startidx) == '-' ...
-                || jsonmodel{i}.rhs(endidx) == '+' ...
-                || jsonmodel{i}.rhs(endidx) == '-'
-            % intercept
-            createdvar = true;
-            if any(strcmp(M_endo_exo_names_trim, 'intercept'))
-                [~, vnames{j}] = fileparts(tempname);
-                vnames{j} = ['intercept_' vnames{j}];
-                assert(~any(strcmp(M_endo_exo_names_trim, vnames{j})));
-            else
-                vnames{j} = 'intercept';
-            end
-            splitstrings{j} = vnames{j};
+    while ~isempty(plus_node) || ~isempty(last_node_to_parse)
+        Xtmp = dseries();
+        if isempty(last_node_to_parse)
+            [plus_node, node_to_parse, last_node_to_parse] = findNextplus_node(plus_node, jsonmodel{i}.line);
         else
-            error('dyn_ols: Shouldn''t arrive here');
+            node_to_parse = last_node_to_parse;
+            last_node_to_parse = [];
         end
-        if createdvar
-            if jsonmodel{i}.rhs(startidx) == '-'
-                Xtmp = dseries(-ones(ds.nobs, 1), ds.firstdate, vnames{j});
+        if strcmp(node_to_parse.node_type, 'VariableNode') || strcmp(node_to_parse.node_type, 'UnaryOpNode')
+            if strcmp(node_to_parse.node_type, 'VariableNode')
+                if strcmp(node_to_parse.type, 'parameter')
+                    % Intercept
+                    Xtmp = dseries(ones(ds.nobs, 1), ds.firstdate, node_to_parse.name);
+                elseif strcmp(node_to_parse.type, 'exogenous') && ~any(strcmp(ds.name, node_to_parse.name))
+                    % Residual if not contained in ds
+                    if isempty(residual)
+                        residual = node_to_parse.name;
+                    else
+                        ols_error(['only one residual allowed per equation; encountered ' residual ' & ' node_to_parse.name], jsonmodel{i}.line);
+                    end
+                elseif strcmp(node_to_parse.type, 'endogenous') ...
+                        || (strcmp(node_to_parse.type, 'exogenous') && any(strcmp(ds.name, node_to_parse.name)))
+                    % Subtract VariableNode from LHS
+                    % NB: treat exogenous that exist in ds as endogenous
+                    lhssub = lhssub + evalNode(ds, node_to_parse, jsonmodel{i}.line, dseries());
+                else
+                    ols_error('unexpected variable type found', jsonmodel{i}.line);
+                end
             else
-                Xtmp = dseries(ones(ds.nobs, 1), ds.firstdate, vnames{j});
+                % Subtract UnaryOpNode from LHS
+                % NB: treat exogenous that exist in ds as endogenous
+                lhssub = lhssub + evalNode(ds, node_to_parse, jsonmodel{i}.line, dseries());
+            end
+        elseif strcmp(node_to_parse.node_type, 'BinaryOpNode') && strcmp(node_to_parse.op, '*')
+            % Parse param_expr * endog_expr
+            Xtmp = parseTimesNode(ds, node_to_parse, jsonmodel{i}.line);
+            if length(Xtmp.name) > 1
+                % Handle constraits
+                % Look through Xtmp names for constant
+                % if found, subtract from LHS
+                to_remove = [];
+                for j = 1:length(Xtmp.name)
+                    if ~isnan(str2double(Xtmp.name{j}))
+                        lhssub = lhssub + str2double(Xtmp.name{j}) * Xtmp{j};
+                        to_remove = [to_remove j];
+                    end
+                    % Multiply by -1 now so that it can be added together below
+                    % Otherwise, it would matter which was encountered first,
+                    % a parameter on its own or a linear constraint
+                    Xtmp.(Xtmp.name{j}) = -1 * Xtmp{j};
+                end
+                for j = length(to_remove):-1:1
+                    Xtmp = Xtmp.remove(Xtmp.name{j});
+                end
             end
-        else
-            Xtmp = eval(regexprep(vnames{j}, regex, 'ds.$&'));
-            Xtmp.rename_(vnames{j});
         end
-        X = [X Xtmp];
+        if ~isempty(Xtmp)
+            to_remove = [];
+            for j = 1:length(Xtmp.name)
+                % Handle constraits
+                idx = find(strcmp(X.name, Xtmp.name{j}));
+                if ~isempty(idx)
+                    X.(X.name{idx}) = X{idx} + Xtmp{j};
+                    to_remove = [to_remove j];
+                end
+            end
+            for j = length(to_remove):-1:1
+                Xtmp = Xtmp.remove(Xtmp.name{j});
+            end
+            X = [X Xtmp];
+        end
     end
 
-    lhssub = getRhsToSubFromLhs(ds, jsonmodel{i}.rhs, regex, splitstrings, pnames);
-    residuals = setdiff(intersect(rhs_, M_.exo_names), ds.name);
-    assert(~isempty(residuals), ['No residuals in equation ' num2str(i)]);
-    assert(length(residuals) == 1, ['More than one residual in equation ' num2str(i)]);
+    clear residual Xtmp plus_node last_node_to_parse to_remove
 
-    Y = eval(regexprep(jsonmodel{i}.lhs, regex, 'ds.$&'));
-    if ~isempty(lhssub)
-        Y = Y - lhssub;
+    Y = Y - lhssub;
+
+    %% Check to see if called from Gibbs
+    st = dbstack(1);
+    varargout = cell(1, 1);
+    called_from_olsgibbs = false;
+    if strcmp(st(1).name, 'olsgibbs')
+        varargout = cell(1, 4);
+        called_from_olsgibbs = true;
     end
 
+    %% Set dates
     fp = max(Y.firstobservedperiod, X.firstobservedperiod);
     lp = min(Y.lastobservedperiod, X.lastobservedperiod);
+    if ~isempty(lhssub)
+        fp = max(fp, lhssub.firstobservedperiod);
+        lp = min(lp, lhssub.lastobservedperiod);
+    end
     if isfield(jsonmodel{i}, 'tags') ...
             && isfield(jsonmodel{i}.tags, 'sample') ...
             && ~isempty(jsonmodel{i}.tags.sample)
@@ -171,7 +206,7 @@ for i = 1:length(jsonmodel)
             fp = fsd;
         end
         if lp < lsd
-             warning(['The sample over which you want to estimate contains NaNs. '...
+            warning(['The sample over which you want to estimate contains NaNs. '...
                 'Adjusting estimation range to end on: ' lp.char])
         else
             lp = lsd;
@@ -179,17 +214,11 @@ for i = 1:length(jsonmodel)
     end
 
     Y = Y(fp:lp);
-    X = X(fp:lp).data;
     if ~isempty(lhssub)
         lhssub = lhssub(fp:lp);
     end
-
-    if isfield(jsonmodel{i}, 'tags') && ...
-            isfield(jsonmodel{i}.tags, 'name')
-        tag = jsonmodel{i}.tags.('name');
-    else
-        tag = ['eq_line_no_' num2str(jsonmodel{i}.line)];
-    end
+    pnames = X.name;
+    X = X(fp:lp).data;
 
     [nobs, nvars] = size(X);
     if called_from_olsgibbs
@@ -203,9 +232,15 @@ for i = 1:length(jsonmodel)
         return
     end
 
+    if isfield(jsonmodel{i}, 'tags') && ...
+            isfield(jsonmodel{i}.tags, 'name')
+        tag = jsonmodel{i}.tags.('name');
+    else
+        tag = ['eq_line_no_' num2str(jsonmodel{i}.line)];
+    end
+
     %% Estimation
     % From LeSage, James P. "Applied Econometrics using MATLAB"
-    
     oo_.ols.(tag).dof = nobs - nvars;
 
     % Estimated Parameters
@@ -214,8 +249,10 @@ for i = 1:length(jsonmodel)
     oo_.ols.(tag).beta = r\(q'*Y.data);
     oo_.ols.(tag).param_idxs = zeros(length(pnames), 1);
     for j = 1:length(pnames)
-        oo_.ols.(tag).param_idxs(j) = find(strcmp(M_.param_names, pnames{j}));
-        M_.params(oo_.ols.(tag).param_idxs(j)) = oo_.ols.(tag).beta(j);
+        if ~strcmp(pnames{j}, 'intercept')
+            oo_.ols.(tag).param_idxs(j) = find(strcmp(M_.param_names, pnames{j}));
+            M_.params(oo_.ols.(tag).param_idxs(j)) = oo_.ols.(tag).beta(j);
+        end
     end
 
     % Yhat
@@ -233,11 +270,9 @@ for i = 1:length(jsonmodel)
     oo_.ols.(tag).resid = Y - oo_.ols.(tag).Yhat;
 
     % Correct Yhat reported back to user
-    if ~isempty(lhssub)
-        Y = Y + lhssub;
-        oo_.ols.(tag).Yhat = oo_.ols.(tag).Yhat + lhssub;
-    end
- 
+    Y = Y + lhssub;
+    oo_.ols.(tag).Yhat = oo_.ols.(tag).Yhat + lhssub;
+
     % Apply correcting function for Yhat if it was passed
     if any(idx) ...
             && length(fitted_names_dict(idx, :)) == 3 ...
@@ -287,12 +322,169 @@ for i = 1:length(jsonmodel)
             sprintf('s^2: %f', oo_.ols.(tag).s2), ...
             sprintf('Durbin-Watson: %f', oo_.ols.(tag).dw)};
 
-        dyn_table(title, preamble, afterward, vnames, ...
-            {'Coefficients','t-statistic','Std. Error'}, 4, ...
+        dyn_table(title, preamble, afterward, pnames, ...
+            {'Estimates','t-statistic','Std. Error'}, 4, ...
             [oo_.ols.(tag).beta oo_.ols.(tag).tstat oo_.ols.(tag).stderr]);
     end
-end
-if ~called_from_olsgibbs
-    varargout{1} = ds;
+
+    if ~called_from_olsgibbs
+        varargout{1} = ds;
+    end
+end
+end
+
+%% Helper Functions
+function ols_error(msg, line)
+error(['ERROR encountered in `dyn_ols` line ' num2str(line) ': ' msg])
+end
+
+function [next_plus_node, node_to_parse, last_node_to_parse] = findNextplus_node(plus_node, line)
+% Given an additive entry in the AST, find the next additive entry
+% (next_plus_node). Also find the node that will be parsed into
+% parameter*endogenous||param||exog|endog (node_to_parse).
+% Function used for moving through the AST.
+if ~(strcmp(plus_node.node_type, 'BinaryOpNode') && strcmp(plus_node.op, '+'))
+    ols_error('pairs of nodes must be separated additively', line);
+end
+next_plus_node = [];
+last_node_to_parse = [];
+if strcmp(plus_node.arg1.node_type, 'BinaryOpNode') && strcmp(plus_node.arg1.op, '+')
+    next_plus_node = plus_node.arg1;
+    node_to_parse = getOlsNode(plus_node.arg2, line);
+elseif strcmp(plus_node.arg2.node_type, 'BinaryOpNode') && strcmp(plus_node.arg2.op, '+')
+    next_plus_node = plus_node.arg2;
+    node_to_parse = getOlsNode(plus_node.arg1, line);
+else
+    node_to_parse = getOlsNode(plus_node.arg1, line);
+    last_node_to_parse = getOlsNode(plus_node.arg2, line);
+end
+end
+
+function node_to_parse = getOlsNode(node, line)
+if ~(strcmp(node.node_type, 'BinaryOpNode') && strcmp(node.op, '*')) ...
+        && ~strcmp(node.node_type, 'VariableNode') ...
+        && ~strcmp(node.node_type, 'UnaryOpNode')
+    ols_error('couldn''t find node to parse', line);
+end
+node_to_parse = node;
+end
+
+function X = parseTimesNode(ds, node, line)
+% Separate the parameter expression from the endogenous expression
+assert(strcmp(node.node_type, 'BinaryOpNode') && strcmp(node.op, '*'))
+[param, X] = parseTimesNodeHelper(ds, node.arg1, line, {}, dseries());
+[param, X] = parseTimesNodeHelper(ds, node.arg2, line, param, X);
+X = X.rename(param{1});
+for ii = 2:length(param)
+    X = [X dseries(X{1}.data, X{1}.firstdate, param{ii})];
+end
+end
+
+function [param, X] = parseTimesNodeHelper(ds, node, line, param, X)
+if isOlsParamExpr(node, line)
+    param = assignParam(param, node, line);
+elseif isOlsVarExpr(node, line)
+    if isempty(X)
+        X = evalNode(ds, node, line, X);
+    else
+        ols_error(['got endog * endog' node.name ' (' node.type ')'], line);
+    end
+else
+    ols_error('unexpected expression', line);
+end
+end
+
+function param = assignParam(param, node, line)
+if ~isempty(param)
+    ols_error(['got param * param' node.name ' (' node.type ')'], line);
+end
+param = assignParamHelper(param, node, line);
+end
+
+function param = assignParamHelper(param, node, line)
+if strcmp(node.node_type, 'NumConstNode')
+    param{end+1} = num2str(node.value);
+elseif strcmp(node.node_type, 'VariableNode')
+    param{end+1} = node.name;
+elseif strcmp(node.node_type, 'BinaryOpNode')
+    if ~strcmp(node.op, '-')
+        ols_error(['got unexpected parameter op ' node.op], line);
+    end
+    param = assignParamHelper(param, node.arg1, line);
+    param = assignParamHelper(param, node.arg2, line);
+else
+    ols_error(['got unexpected node (' node.type ')'], line);
+end
+end
+
+function tf = isOlsVar(node)
+if strcmp(node.node_type, 'VariableNode') ...
+        && (strcmp(node.type, 'endogenous') ...
+        || (strcmp(node.type, 'exogenous') && any(strcmp(ds.name, node.name))))
+    tf = true;
+elseif strcmp(node.node_type, 'UnaryOpNode') ...
+        && (strcmp(node.arg.type, 'endogenous') ...
+        || (strcmp(node.arg.type, 'exogenous') && any(strcmp(ds.name, node.arg.name))))
+    tf = true;
+else
+    tf = false;
+end
+end
+
+function tf = isOlsVarExpr(node, line)
+if strcmp(node.node_type, 'VariableNode') || strcmp(node.node_type, 'UnaryOpNode')
+    tf = isOlsVar(node);
+elseif strcmp(node.node_type, 'BinaryOpNode')
+    tf = isOlsVarExpr(node.arg1, line) || isOlsVarExpr(node.arg2, line);
+else
+    ols_error(['got unexpected type ' node.node_type], line);
+end
+end
+
+function X = evalNode(ds, node, line, X)
+if strcmp(node.node_type, 'NumConstNode')
+    X = dseries(str2double(node.value), ds.dates, 'const');
+elseif strcmp(node.node_type, 'VariableNode')
+    if ~(strcmp(node.type, 'endogenous') ...
+            || (strcmp(node.type, 'exogenous') && any(strcmp(ds.name, node.name))))
+        ols_error(['got unexpected type ' node.name ': ' node.type], line);
+    end
+    X = ds.(node.name)(node.lag);
+elseif strcmp(node.node_type, 'UnaryOpNode')
+    Xtmp = evalNode(ds, node.arg, line, X);
+    % Only works if dseries supports . notation for unary op (true for log/diff)
+    % Otherwise, use: X = eval([node.op '(Xtmp)']);
+    X = Xtmp.(node.op);
+elseif strcmp(node.node_type, 'BinaryOpNode')
+    Xtmp1 = evalNode(ds, node.arg1, line, X);
+    Xtmp2 = evalNode(ds, node.arg2, line, X);
+    X = X + eval(['Xtmp1 ' node.op ' Xtmp2']);
+else
+    ols_error(['got unexpected node type ' node.node_type], line);
+end
+end
+
+function tf = isOlsParam(node)
+if strcmp(node.node_type, 'VariableNode') && strcmp(node.type, 'parameter')
+    tf = true;
+else
+    tf = false;
+end
+end
+
+function tf = isOlsParamExpr(node, line)
+if strcmp(node.node_type, 'NumConstNode')
+    tf = true;
+elseif strcmp(node.node_type, 'VariableNode')
+    tf = isOlsParam(node);
+elseif strcmp(node.node_type, 'UnaryOpNode')
+    tf = false;
+elseif strcmp(node.node_type, 'BinaryOpNode')
+    tf = isOlsParamExpr(node.arg1) && isOlsParamExpr(node.arg2);
+    if tf && ~strcmp(node.op, '-')
+        ols_error(['got unexpected op ' node.op], line);
+    end
+else
+    ols_error(['got unexpected type ' node.node_type], line);
 end
 end
diff --git a/tests/estimation/univariate/ols/ols.mod b/tests/estimation/univariate/ols/ols.mod
new file mode 100644
index 000000000..ca4b037ee
--- /dev/null
+++ b/tests/estimation/univariate/ols/ols.mod
@@ -0,0 +1,150 @@
+// --+ options: json=compute nopreprocessoroutput stochastic +--
+
+/* REMARK
+** ------
+**
+** You need to have the first line on top of the mod file. The options defined on this line are passed
+** to the dynare command (you can add other options, separated by spaces or commas). The option defined
+** here is mandatory for the decomposition. It forces Dynare to output another representation of the
+** model in JSON file (additionaly to the matlab files) which is used here to manipulate the equations.
+*/
+
+var
+U2_Q_YED
+U2_G_YER
+U2_STN
+U2_ESTN
+U2_EHIC
+DE_Q_YED
+DE_G_YER
+DE_EHIC
+
+;
+
+varexo
+res_U2_Q_YED
+res_U2_G_YER
+res_U2_STN
+res_U2_ESTN
+res_U2_EHIC
+res_DE_Q_YED
+res_DE_G_YER
+res_DE_EHIC
+;
+
+parameters
+u2_q_yed_ecm_u2_q_yed_L1
+u2_q_yed_ecm_u2_stn_L1
+u2_q_yed_u2_g_yer_L1
+u2_q_yed_u2_stn_L1
+u2_g_yer_ecm_u2_q_yed_L1
+u2_g_yer_ecm_u2_stn_L1
+u2_g_yer_u2_q_yed_L1
+u2_g_yer_u2_g_yer_L1
+u2_g_yer_u2_stn_L1
+u2_stn_ecm_u2_q_yed_L1
+u2_stn_ecm_u2_stn_L1
+u2_stn_u2_q_yed_L1
+u2_stn_u2_g_yer_L1
+u2_estn_u2_estn_L1
+u2_ehic_u2_ehic_L1
+
+de_q_yed_ecm_de_q_yed_L1
+de_q_yed_ecm_u2_stn_L1
+de_q_yed_de_g_yer_L1
+de_q_yed_u2_stn_L1
+de_g_yer_ecm_de_q_yed_L1
+de_g_yer_ecm_u2_stn_L1
+de_g_yer_de_q_yed_L1
+de_g_yer_de_g_yer_L1
+de_g_yer_u2_stn_L1
+de_ehic_de_ehic_L1
+
+
+;
+
+u2_q_yed_ecm_u2_q_yed_L1  = -0.82237516589315   ;
+u2_q_yed_ecm_u2_stn_L1    = -0.323715338568976  ;
+u2_q_yed_u2_g_yer_L1      =  0.0401361895021084 ;
+u2_q_yed_u2_stn_L1        =  0.058397703958446  ;
+u2_g_yer_ecm_u2_q_yed_L1  =  0.0189896046977421 ;
+u2_g_yer_ecm_u2_stn_L1    = -0.109597659887432  ;
+u2_g_yer_u2_q_yed_L1      =  0.0037667967632025 ;
+u2_g_yer_u2_g_yer_L1      =  0.480506381923644  ;
+u2_g_yer_u2_stn_L1        = -0.0722359286123494 ;
+u2_stn_ecm_u2_q_yed_L1    = -0.0438500662608356 ;
+u2_stn_ecm_u2_stn_L1      = -0.153283917138772  ;
+u2_stn_u2_q_yed_L1        =  0.0328744983772825 ;
+u2_stn_u2_g_yer_L1        =  0.292121949736756  ;
+u2_estn_u2_estn_L1        =  1                  ;
+u2_ehic_u2_ehic_L1        =  1                  ;
+
+de_q_yed_ecm_de_q_yed_L1  = -0.822375165893149  ;
+de_q_yed_ecm_u2_stn_L1    = -0.323715338568977  ;
+de_q_yed_de_g_yer_L1      =  0.0401361895021082 ;
+de_q_yed_u2_stn_L1        =  0.0583977039584461 ;
+de_g_yer_ecm_de_q_yed_L1  =  0.0189896046977422 ;
+de_g_yer_ecm_u2_stn_L1    = -0.109597659887433  ;
+de_g_yer_de_q_yed_L1      =  0.00376679676320256;
+de_g_yer_de_g_yer_L1      =  0.480506381923643  ;
+de_g_yer_u2_stn_L1        = -0.0722359286123494 ;
+de_ehic_de_ehic_L1        =  1                  ;
+
+
+model(linear);
+[name = 'eq1']
+diff(U2_Q_YED) =   u2_q_yed_ecm_u2_q_yed_L1 * (U2_Q_YED(-1) - U2_EHIC(-1))
+                 + u2_q_yed_ecm_u2_stn_L1   * (U2_STN(-1)   - U2_ESTN(-1))
+                 + u2_q_yed_u2_g_yer_L1     * diff(U2_G_YER(-1))
+                 + u2_q_yed_u2_stn_L1       * diff(U2_STN(-1))
+                 + res_U2_Q_YED                                           ;
+[name = 'eq2']
+diff(U2_G_YER) =   u2_g_yer_ecm_u2_q_yed_L1 * (U2_Q_YED(-1) - U2_EHIC(-1))
+                 + u2_g_yer_ecm_u2_stn_L1   * (U2_STN(-1)   - U2_ESTN(-1))
+                 + u2_g_yer_u2_q_yed_L1     * diff(U2_Q_YED(-1))
+                 + u2_g_yer_u2_g_yer_L1     * diff(U2_G_YER(-1))
+                 + u2_g_yer_u2_stn_L1       * diff(U2_STN(-1))
+                 + res_U2_G_YER                                           ;
+[name = 'eq3']
+diff(U2_STN)   =   u2_stn_ecm_u2_q_yed_L1   * (U2_Q_YED(-1) - U2_EHIC(-1))
+                 + u2_stn_ecm_u2_stn_L1     * (U2_STN(-1)   - U2_ESTN(-1))
+                 + u2_stn_u2_q_yed_L1       * diff(U2_Q_YED(-1))
+                 + u2_stn_u2_g_yer_L1       * diff(U2_G_YER(-1))
+                 + res_U2_STN                                             ;
+[name = 'eq4']
+U2_ESTN        =   u2_estn_u2_estn_L1       * U2_ESTN(-1)
+                 + res_U2_ESTN                                            ;
+[name = 'eq5']
+U2_EHIC        =   u2_ehic_u2_ehic_L1       * U2_EHIC(-1)
+                 + res_U2_EHIC                                            ;
+[name = 'eq6']
+diff(DE_Q_YED) =   de_q_yed_ecm_de_q_yed_L1 * (DE_Q_YED(-1) - DE_EHIC(-1))
+                 + de_q_yed_ecm_u2_stn_L1   * (U2_STN(-1)   - U2_ESTN(-1))
+                 + de_q_yed_de_g_yer_L1     * diff(DE_G_YER(-1))
+                 + de_q_yed_u2_stn_L1       * diff(U2_STN(-1))
+                 + res_DE_Q_YED                                           ;
+[name = 'eq7']
+diff(DE_G_YER) =   de_g_yer_ecm_de_q_yed_L1 * (DE_Q_YED(-1) - DE_EHIC(-1))
+                 + de_g_yer_ecm_u2_stn_L1   * (U2_STN(-1)   - U2_ESTN(-1))
+                 + de_g_yer_de_q_yed_L1     * diff(DE_Q_YED(-1))
+                 + de_g_yer_de_g_yer_L1     * diff(DE_G_YER(-1))
+                 + de_g_yer_u2_stn_L1       * diff(U2_STN(-1))
+                 + res_DE_G_YER                                           ;
+[name = 'eq8']
+DE_EHIC        =   de_ehic_de_ehic_L1       * DE_EHIC(-1)
+                 + res_DE_EHIC                                            ;
+
+
+
+end;
+
+shocks;
+var res_U2_Q_YED = 0.005;
+var res_U2_G_YER = 0.005;
+var res_U2_STN = 0.005;
+var res_U2_ESTN = 0.005;
+var res_U2_EHIC = 0.005;
+var res_DE_Q_YED = 0.005;
+var res_DE_G_YER = 0.005;
+var res_DE_EHIC = 0.005;
+end;
diff --git a/tests/estimation/univariate/ols/run_simulation_test.m b/tests/estimation/univariate/ols/run_simulation_test.m
new file mode 100644
index 000000000..090832c81
--- /dev/null
+++ b/tests/estimation/univariate/ols/run_simulation_test.m
@@ -0,0 +1,39 @@
+close all
+
+dynare ols.mod;
+
+global options_
+options_.noprint = true;
+
+NSIMS = 1000;
+
+calibrated_values = M_.params;
+Sigma_e = M_.Sigma_e;
+
+options_.bnlms.set_dynare_seed_to_default = false;
+
+nparampool = length(M_.params);
+BETA = zeros(NSIMS, nparampool);
+for i=1:NSIMS
+    i
+    firstobs = rand(3, length(M_.endo_names));
+    M_.params = calibrated_values;
+    M_.Sigma_e = Sigma_e; 
+    simdata = simul_backward_model(dseries(firstobs, dates('1995Q1'), M_.endo_names), 10000);
+    simdata = simdata(simdata.dates(5001:6000));
+    names = regexp(simdata.name, 'res\w*');
+    idxs = find(cellfun(@isempty, names));
+    dyn_ols(simdata{idxs}, {}, {'eq7'});
+    BETA(i, :) = M_.params';
+end
+
+mean(BETA)' - calibrated_values
+
+for i=1:nparampool
+    figure
+    hold on
+    title(strrep(M_.param_names(i,:), '_', '\_'));
+    histogram(BETA(:,i),50);
+    line([calibrated_values(i) calibrated_values(i)], [0 NSIMS/10], 'LineWidth', 2, 'Color', 'r');
+    hold off
+end
\ No newline at end of file