diff --git a/tests/ep/burnside.mod b/tests/ep/burnside.mod
new file mode 100644
index 000000000..f6f1fb14d
--- /dev/null
+++ b/tests/ep/burnside.mod
@@ -0,0 +1,104 @@
+addpath ..
+var y x;
+varexo e;
+
+parameters beta theta rho xbar;
+xbar = 0.0179;
+rho =  -0.139;
+theta = -1.5;
+beta = 0.95;
+
+model(use_dll);
+1 = beta*exp(theta*x(+1))*(1+y(+1))/y;
+x = (1-rho)*xbar + rho*x(-1)+e;
+end;
+
+shocks;
+var e; stderr 0.0348;
+end;
+
+initval;
+x = xbar;
+y = beta*exp(theta*xbar)/(1-beta*exp(theta*xbar));
+end;
+
+resid(1);
+
+steady;
+
+if beta*exp(theta*xbar+.5*theta^2*M_.Sigma_e/(1-rho)^2)>1-eps
+   disp('The model doesn''t have a solution!')
+   return
+end
+
+set_dynare_seed('default');
+stoch_simul(order=1,irf=0,periods=5000);
+y_perturbation_1 = oo_.endo_simul(1,:)';
+
+set_dynare_seed('default');
+stoch_simul(order=2,irf=0,periods=5000);
+y_perturbation_2 = oo_.endo_simul(1,:)';
+
+set_dynare_seed('default');
+stoch_simul(order=2,pruning,irf=0,periods=5000);
+y_perturbation_2_pruning = oo_.endo_simul(1,:)';
+
+options_.maxit_ = 100;
+options_.ep.verbosity = 0;
+options_.ep.stochastic.order = 0;
+options_.ep.stochastic.nodes = 2;
+options_.console_mode = 0;
+
+set_dynare_seed('default');
+ts = extended_path([],5000);
+
+options_.ep.stochastic.order = 2;
+set_dynare_seed('default');
+ts1_4 = extended_path([],5000);
+
+set_dynare_seed('default');
+ytrue=exact_solution(M_,oo_,800);
+
+disp('True mean and standard deviation')
+disp(mean(ytrue(101:end)))
+disp(sqrt(var(ytrue(101:end))))
+
+disp('Order 1 perturbation mean and standard deviation')
+disp(mean(y_perturbation_1(101:end)))
+disp(sqrt(var(y_perturbation_1(101:end))))
+
+disp('Order 2 perturbation mean and standard deviation')
+disp(mean(y_perturbation_2(101:end)))
+disp(sqrt(var(y_perturbation_2(101:end))))
+
+disp('Order 2 perturbation (with pruning) mean and standard deviation')
+disp(mean(y_perturbation_2_pruning(101:end)))
+disp(sqrt(var(y_perturbation_2_pruning(101:end))))
+
+disp('Extended path mean and standard deviation')
+disp(mean(ts(1,101:end)))
+disp(sqrt(var(ts(1,101:end))))
+
+disp('Stochastic extended path mean and standard deviation')
+disp(mean(ts1_4(1,101:end)))
+disp(sqrt(var(ts1_4(1,101:end))))
+
+disp('Accuracy error (order 1 perturbation)')
+disp(mean(100*abs(y_perturbation_1-ytrue')./ytrue'));
+disp(max(100*abs(y_perturbation_1-ytrue')./ytrue'));
+
+disp('Accuracy error (order 2 perturbation)')
+disp(mean(100*abs(y_perturbation_2-ytrue')./ytrue'));
+disp(max(100*abs(y_perturbation_2-ytrue')./ytrue'));
+
+disp('Accuracy error (order 2 perturbation with pruning)')
+disp(100*mean(abs(y_perturbation_2_pruning-ytrue')./ytrue'));
+disp(100*max(abs(y_perturbation_2_pruning-ytrue')./ytrue'));
+
+disp('Accuracy error (extended path)')
+disp(mean(100*abs(ts(1,:)'-ytrue')./ytrue'));
+disp(max(100*abs(ts(1,:)'-ytrue')./ytrue'));
+
+disp('Accuracy error (stochastic extended path)')
+disp(mean(100*abs(ts1_4(1,:)'-ytrue')./ytrue'));
+disp(max(100*abs(ts1_4(1,:)'-ytrue')./ytrue'));
diff --git a/tests/ep/exact_solution.m b/tests/ep/exact_solution.m
new file mode 100644
index 000000000..5b525797f
--- /dev/null
+++ b/tests/ep/exact_solution.m
@@ -0,0 +1,30 @@
+function y=exact_solution(M,oo,n)
+    beta = M.params(1);
+    theta = M.params(2);
+    rho = M.params(3);
+    xbar = M.params(4);
+    sigma2 = M.Sigma_e;
+    
+    if beta*exp(theta*xbar+.5*theta^2*sigma2/(1-rho)^2)>1-eps
+        disp('The model doesn''t have a solution!')
+        return
+    end
+    
+    i = 1:n;
+    a = theta*xbar*i+(theta^2*sigma2)/(2*(1-rho)^2)*(i-2*rho*(1-rho.^i)/(1-rho)+rho^2*(1-rho.^(2*i))/(1-rho^2));
+    b = theta*rho*(1-rho.^i)/(1-rho);
+    
+    x = oo.endo_simul(2,:);
+    xhat = x-xbar;
+    
+    n2 = size(x,2);
+    
+    y = zeros(1,n2);
+    
+    
+    for j=1:n2
+        y(j) = sum(beta.^i.*exp(a+b*xhat(j)));
+    end
+    
+    disp(sum(beta.^i.*exp(theta*xbar*i)))
+    disp(sum(beta.^i.*exp(a)))
\ No newline at end of file