function [s, penultimatepoint1, penultimatepoint2,... 
             switchpoint1, switchpoint2,...
             history1, history2] = ...
    mymkdata1(M_,oo_,nperiods,order,active_shocks_labels,active_shocks,...
              init1,init2,output_var_labels)

dr_ = oo_.dr;

%  dr_ structure returned by dynare with model's decision rule
%  nperiods number of periods for simulated data
%  order of the decision rule to be used.


nvars = size(dr_.ghx,1);
nactive_shocks = size(active_shocks_labels,1);
nstatic = M_.nstatic;
nfwrd = M_.nfwrd;
statevar_pos = (nstatic +1):(nvars-nfwrd);
%statevar_pos = (dr_.nstatic +1):(dr_.nstatic+size(dr_.ghx,2));

nshocks = M_.exo_nbr;
active_shock_pos_vec = zeros(nactive_shocks,1);
for this_shock = 1:nactive_shocks
    active_shock_pos_vec(this_shock) = strmatch(active_shocks_labels(this_shock,:),M_.exo_names,'exact');
end

shocks = zeros(nshocks,nperiods);

if ~exist('active_shocks','var')
s1 = RandStream('mt19937ar','seed',1);
shocks(active_shock_pos_vec,:) = randn(s1,nactive_shocks,nperiods);
else 
shocks(active_shock_pos_vec,:) = active_shocks;

end

if ~exist('init1','var')
   init1 = zeros(nvars,1);
   init2 = zeros(nvars,1);
end

history = zeros(nvars,nperiods+1);
history(:,1) = init1(oo_.dr.order_var);

for i = 2:nperiods+1   
    history(:,i) = dr_.ghx*history(statevar_pos,i-1)+ dr_.ghu*shocks(:,i-1);
end
penultimatepoint1 = history(dr_.inv_order_var,end-1);
switchpoint1 = history(dr_.inv_order_var,end);

if order>1
   history2 = zeros(nvars,nperiods+1);
   history2(:,1) = init2(oo_.dr.order_var);
     for i = 2:nperiods+1
     history2(:,i) = 0.5*dr_.ghs2 + dr_.ghx*history2(statevar_pos,i-1) + dr_.ghu*shocks(:,i-1) + ... 
                     0.5*dr_.ghxx*kron(history(statevar_pos,i-1),history(statevar_pos,i-1)) + ...
                     0.5*dr_.ghuu*kron(shocks(:,i-1),shocks(:,i-1)) + ...
                     dr_.ghxu*kron(history(statevar_pos,i-1),shocks(:,i-1));
     end
     
     penultimatepoint2 = history2(dr_.inv_order_var,end-1);
     switchpoint2 = history2(dr_.inv_order_var,end);
     
     %history currently contains the first order responses; reassign it to
     %contain the second-order responses, but also prepare history1 and
     %history2 for output -- rearrange them in order of definition in the
     %.mod file
     history1 = history(dr_.inv_order_var,:);
     history = history2;
     history2 = history2(dr_.inv_order_var,:);
else 
     history1 = history(dr_.inv_order_var,:);
     history2 = nan*history;
     penultimatepoint2 = nan;
     switchpoint2 = nan;
end
 

                                          
% this loop assigns each row of history to its variable
% it also assigns steady state values
% IRFs are stored under varname_irf  where varname is the variable name
% SS values are stored under varname
% NB: SS values are scalar, IRFs are column vectors.


reordered = M_.endo_names(dr_.order_var,:);          % the vector reordered contains names of variables
                                            % corresponding to rows of
                                            % history

history = history';

if ~exist('output_var_labels','var')
    for indxi = 1:M_.endo_nbr
        eval(['s','.',deblank(M_.endo_names(indxi,:)),'_ss = oo_.dr.ys(indxi);']);
        eval(['s','.',deblank(reordered(indxi,:)),'_irf=history(:,indxi);']);
    end
else
    noutput_var_labels = size(output_var_labels,1);
    for indxi=1:noutput_var_labels
        this_var_label = deblank(output_var_labels(indxi,:));
        ss_pos = strmatch(this_var_label,M_.endo_names,'exact');
        eval(['s','.',this_var_label,'_ss = oo_.dr.ys(ss_pos);']);
        
        var_pos = strmatch(this_var_label,reordered,'exact');
        eval(['s','.',this_var_label,'_irf=history(:,var_pos);']);
    end
end

