diff --git a/meson.build b/meson.build
index f8d230fdf8d68bad83f4df8e939b5db488fe2958..19ad92b4c49e45b435010dbb45547c51c95139eb 100644
--- a/meson.build
+++ b/meson.build
@@ -804,10 +804,12 @@ mod_and_m_tests = [
                 'estimation/fsdat_simul.m' ] },
   { 'test' : [ 'estimation/fs2000.mod' ],
     'extra' : [ 'estimation/fsdat_simul.m' ] },
-  { 'test' : [ 'gsa/ls2003a.mod' ],
-    'extra' : [ 'gsa/data_ca1.m' ] },
-  { 'test' : [ 'gsa/ls2003.mod' ],
-    'extra' : [ 'gsa/data_ca1.m' ] },
+  { 'test' : [ 'gsa/ls2003a.mod',
+               'gsa/ls2003.mod',
+               'gsa/ls2003scr.mod',
+               'gsa/ls2003ide.mod' ],
+    'extra' : [ 'gsa/data_ca1.m',
+                'gsa/ls2003_model.inc' ] },
   { 'test' : [ 'gsa/cod_ML_morris/cod_ML_morris.mod' ] },
   { 'test' : [ 'gsa/morris/nk_est.mod' ],
     'extra' : [ 'gsa/morris/nk_est_data.m' ] },
diff --git a/tests/gsa/ls2003.mod b/tests/gsa/ls2003.mod
index 25e02deffe02f9b97a50305bc10b452d59b3682c..9443b535f1b8154ef7f2c15bf745d9f299c09061 100644
--- a/tests/gsa/ls2003.mod
+++ b/tests/gsa/ls2003.mod
@@ -1,45 +1,4 @@
-var y y_s R pie $\pi$ dq pie_s de A y_obs pie_obs R_obs ${R^{obs}}$;
-varexo e_R ${e^R}$ e_q e_ys e_pies e_A;
-
-parameters psi1 psi2 psi3 rho_R ${\rho_R}$  tau alpha rr k rho_q rho_A rho_ys rho_pies;
-
-psi1 = 1.54;
-psi2 = 0.25;
-psi3 = 0.25;
-rho_R = 0.5;
-alpha = 0.3;
-rr = 2.51;
-k = 0.5;
-tau = 0.5;
-rho_q = 0.4;
-rho_A = 0.2;
-rho_ys = 0.9;
-rho_pies = 0.7;
-
-
-model(linear);
-y = y(+1) - (tau +alpha*(2-alpha)*(1-tau))*(R-pie(+1))-alpha*(tau +alpha*(2-alpha)*(1-tau))*dq(+1) + alpha*(2-alpha)*((1-tau)/tau)*(y_s-y_s(+1))-A(+1);
-pie = exp(-rr/400)*pie(+1)+alpha*exp(-rr/400)*dq(+1)-alpha*dq+(k/(tau+alpha*(2-alpha)*(1-tau)))*y+k*alpha*(2-alpha)*(1-tau)/(tau*(tau+alpha*(2-alpha)*(1-tau)))*y_s;
-pie = de+(1-alpha)*dq+pie_s;
-R = rho_R*R(-1)+(1-rho_R)*(psi1*pie+psi2*(y+alpha*(2-alpha)*((1-tau)/tau)*y_s)+psi3*de)+e_R;
-dq = rho_q*dq(-1)+e_q;
-y_s = rho_ys*y_s(-1)+e_ys;
-pie_s = rho_pies*pie_s(-1)+e_pies;
-A = rho_A*A(-1)+e_A;
-y_obs = y-y(-1)+A;
-pie_obs = 4*pie;
-R_obs = 4*R;
-end;
-
-shocks;
-	var e_R = 1.25^2;
-	var e_q = 2.5^2;
-	var e_A = 1.89;
-	var e_ys = 1.89;
-	var e_pies = 1.89;
-end;
-
-varobs y_obs R_obs pie_obs dq de;
+@#include "ls2003_model.inc"
 
 estimated_params;
 psi1 , gamma_pdf,1.5,0.5;
@@ -66,7 +25,7 @@ options_.TeX=true;
 disp(' ');
 disp('NOW I DO STABILITY MAPPING and prepare sample for Reduced form Mapping');
 disp(' ');
-disp('Press ENTER to continue'); pause(5);
+% disp('Press ENTER to continue'); pause(5);
 
 dynare_sensitivity(redform=1,nodisplay,Nsam=512); //create sample of reduced form coefficients
 // NOTE: since namendo is empty by default, 
@@ -75,7 +34,7 @@ dynare_sensitivity(redform=1,nodisplay,Nsam=512); //create sample of reduced for
 disp(' ');
 disp('ANALYSIS OF REDUCED FORM COEFFICIENTS');
 disp(' ');
-disp('Press ENTER to continue'); pause(5);
+% disp('Press ENTER to continue'); pause(5);
 
 dynare_sensitivity(nodisplay, load_stab=1,  // load previously generated sample analysed for stability
 redform=1,  // do the reduced form mapping
@@ -87,6 +46,7 @@ stab=0, // don't repeat again the stability mapping
 Nsam=512);
 
 
+
 disp(' ');
 disp('THE PREVIOUS TWO CALLS COULD BE DONE TOGETHER');
 disp('BY USING THE COMBINED CALL');
@@ -94,7 +54,7 @@ disp(' ');
 disp('dynare_sensitivity(redform=1,')
 disp('threshold_redform=[-1 0], namendo=(pie,R), namexo=(e_R), namlagendo=(R));')
 disp(' ');
-disp('Press ENTER to continue'); pause(5);
+% disp('Press ENTER to continue'); pause(5);
 //dynare_sensitivity(
 //redform=1, //create sample of reduced form coefficients
 //logtrans_redform=1,  // estimate log-transformed reduced form coefficients (default=0)
@@ -107,7 +67,7 @@ disp('Press ENTER to continue'); pause(5);
 
 disp(' ');
 disp('MC FILTERING(rmse=1), TO MAP THE FIT FROM PRIORS');
-disp('Press ENTER to continue'); pause(5);
+%disp('Press ENTER to continue'); pause(5);
 
 dynare_sensitivity(nodisplay, datafile='data_ca1.m',first_obs=8,nobs=79,prefilter=1, // also presample=2,loglinear, are admissible
 load_stab=1,     // load prior sample
@@ -125,7 +85,7 @@ disp('logtrans_redform=1, namendo=(pie,R), namexo=(e_R), namlagendo=(R),')
 disp('datafile=data_ca1.m,first_obs=8,nobs=79,prefilter=1,')
 disp('istart_rmse=2, rmse=1);')
 disp(' ');
-disp('Press ENTER to continue'); pause(5);
+% disp('Press ENTER to continue'); pause(5);
 //dynare_sensitivity(
 //redform=1, //create sample of reduced form coefficients
 //logtrans_redform=1,  // estimate log-transformed reduced form coefficients (default=0)
@@ -142,14 +102,13 @@ disp('Press ENTER to continue'); pause(5);
 disp(' ');
 disp('I ESTIMATE THE MODEL');
 disp(' ');
-disp('Press ENTER to continue'); pause(5);
+% disp('Press ENTER to continue'); pause(5);
 
 // run this to generate posterior mode and Metropolis files if not yet done
 estimation(datafile='data_ca1.m',first_obs=8,nobs=79,mh_nblocks=1,
    prefilter=1,mh_jscale=0.5,mh_replic=5000,silent_optimizer, mode_compute=4, mh_drop=0.6, nodisplay,
    bayesian_irf, filtered_vars, smoother) y_obs R_obs pie_obs dq de;
 
-
 // run this to produce posterior samples of filtered, smoothed and irf variables, if not yet done
 //estimation(datafile='data_ca1.m',first_obs=8,nobs=79,mh_nblocks=2,prefilter=1,mh_jscale=0.3,
 //          mh_replic=0, mode_file='ls2003/Output/ls2003_mode', mode_compute=0, load_mh_file, bayesian_irf,
@@ -160,7 +119,7 @@ disp('WE DO STABILITY MAPPING AGAIN, EXPLORING AROUND THE MODE BY A WIDTH +- 0.2
 disp('THEN WE DO RMSE analysis in this neighborhood');
 disp('Typical for ML estimation, also feasible for posterior mode');
 disp(' ');
-disp('Press ENTER to continue'); pause(5);
+% disp('Press ENTER to continue'); pause(5);
 
 dynare_sensitivity(nodisplay, pprior=0,Nsam=512,neighborhood_width=0.2,
 mode_file='ls2003/Output/ls2003_mode',  // specifies the mode file where the mode and Hessian are stored
@@ -171,7 +130,7 @@ disp(' ');
 disp('WE DO STABILITY MAPPING AGAIN, BUT FOR MULTIVARIATE SAMPLE AT THE POSTERIOR MODE (or ML) and Hessian (pprior=0 & ppost=0)');
 disp('Typical for ML estimation, also feasible for posterior mode');
 disp(' ');
-disp('Press ENTER to continue'); pause(5);
+% disp('Press ENTER to continue'); pause(5);
 
 dynare_sensitivity(nodisplay, pprior=0,Nsam=512,
 mode_file='ls2003/Output/ls2003_mode'  // specifies the mode file where the mode and Hessian are stored
@@ -181,7 +140,7 @@ mode_file='ls2003/Output/ls2003_mode'  // specifies the mode file where the mode
 disp(' ');
 disp('RMSE ANALYSIS FOR MULTIVARIATE SAMPLE AT THE POSTERIOR MODE');
 disp(' ');
-disp('Press ENTER to continue'); pause(5);
+disp('Press ENTER to continue'); % pause(5);
 dynare_sensitivity(nodisplay, mode_file='ls2003/Output/ls2003_mode',
 datafile='data_ca1.m',first_obs=8,nobs=79,prefilter=1,
 pprior=0,
@@ -199,10 +158,10 @@ disp('dynare_sensitivity(pprior=0,Nsam=512,alpha2_stab=0.4,mode_file=ls2003_mode
 disp('datafile=data_ca1.m,first_obs=8,nobs=79,prefilter=1,')
 disp('rmse=1, alpha2_rmse=0, alpha_rmse=0);')
 disp(' ');
-disp('Press ENTER to continue'); pause(5);
+% disp('Press ENTER to continue'); pause(5);
 //dynare_sensitivity(pprior=0,Nsam=2048,alpha2_stab=0.4,mode_file='ls2003/Output/ls2003_mode',
 //datafile='data_ca1.m',first_obs=8,nobs=79,prefilter=1,
-//rmse=1,
+//rmse=1pae
 //alpha2_rmse=0, // no correlation analysis
 //alpha_rmse=0  // no Smirnov sensitivity analysis
 //);
@@ -210,7 +169,7 @@ disp('Press ENTER to continue'); pause(5);
 disp(' ');
 disp('RMSE ANALYSIS FOR POSTERIOR MCMC sample (ppost=1)');
 disp('Needs a call to dynare_estimation to load all MH environment');
-disp('Press ENTER to continue'); pause(5);
+% disp('Press ENTER to continue'); pause(5);
 //estimation(datafile='data_ca1.m',first_obs=8,nobs=79,mh_nblocks=2, mode_file='ls2003/Output/ls2003_mode', load_mh_file,
 //  prefilter=1,mh_jscale=0.5,mh_replic=0, mode_compute=0, mh_drop=0.6);
 
diff --git a/tests/gsa/ls2003_model.inc b/tests/gsa/ls2003_model.inc
new file mode 100644
index 0000000000000000000000000000000000000000..c8e3fea8a53fbbfb35264ef8890a47af8a696383
--- /dev/null
+++ b/tests/gsa/ls2003_model.inc
@@ -0,0 +1,70 @@
+var y 
+    y_s 
+    R 
+    pie $\pi$ 
+    dq  
+    pie_s 
+    de 
+    A   
+    y_obs   ${y^{obs}}$ 
+    pie_obs ${\pi^{obs}}$
+    R_obs ${R^{obs}}$;
+
+varexo e_R ${\varepsilon^R}$ 
+    e_q ${\varepsilon^q}$ 
+    e_ys ${\varepsilon^{ys}}$
+    e_pies ${\varepsilon^\pi}$
+    e_A ${\varepsilon^A}$;
+
+parameters psi1 ${\psi_1}$ 
+    psi2 ${\psi_2}$ 
+    psi3 ${\psi_3}$
+    rho_R ${\rho_R}$  
+    tau    ${\tau}$
+    alpha ${\alpha}$
+    rr 
+    k 
+    rho_q ${\rho_q}$
+    rho_A ${\rho_A}$
+    rho_ys ${\rho_{ys}}$
+    rho_pies    ${\rho_\pi}$;
+
+psi1 = 1.54;
+psi2 = 0.25;
+psi3 = 0.25;
+rho_R = 0.5;
+alpha = 0.3;
+rr = 2.51;
+k = 0.5;
+tau = 0.5;
+rho_q = 0.4;
+rho_A = 0.2;
+rho_ys = 0.9;
+rho_pies = 0.7;
+
+
+model(linear);
+y = y(+1) - (tau +alpha*(2-alpha)*(1-tau))*(R-pie(+1))-alpha*(tau +alpha*(2-alpha)*(1-tau))*dq(+1) + alpha*(2-alpha)*((1-tau)/tau)*(y_s-y_s(+1))-A(+1);
+pie = exp(-rr/400)*pie(+1)+alpha*exp(-rr/400)*dq(+1)-alpha*dq+(k/(tau+alpha*(2-alpha)*(1-tau)))*y+k*alpha*(2-alpha)*(1-tau)/(tau*(tau+alpha*(2-alpha)*(1-tau)))*y_s;
+pie = de+(1-alpha)*dq+pie_s;
+R = rho_R*R(-1)+(1-rho_R)*(psi1*pie+psi2*(y+alpha*(2-alpha)*((1-tau)/tau)*y_s)+psi3*de)+e_R;
+dq = rho_q*dq(-1)+e_q;
+y_s = rho_ys*y_s(-1)+e_ys;
+pie_s = rho_pies*pie_s(-1)+e_pies;
+A = rho_A*A(-1)+e_A;
+y_obs = y-y(-1)+A;
+pie_obs = 4*pie;
+R_obs = 4*R;
+end;
+
+shocks;
+	var e_R = 1.25^2;
+	var e_q = 2.5^2;
+	var e_A = 1.89;
+	var e_ys = 1.89;
+	var e_pies = 1.89;
+end;
+
+varobs y_obs R_obs pie_obs dq de;
+
+options_.TeX=1;
\ No newline at end of file
diff --git a/tests/gsa/ls2003a.mod b/tests/gsa/ls2003a.mod
index 355d332d1e24a9181df9835ac5aea1162e6c9949..a25e062e8ff843c9c9188f255948b187efe3722e 100644
--- a/tests/gsa/ls2003a.mod
+++ b/tests/gsa/ls2003a.mod
@@ -1,45 +1,4 @@
-var y y_s R pie dq pie_s de A y_obs pie_obs R_obs;
-varexo e_R e_q e_ys e_pies e_A;
-
-parameters psi1 psi2 psi3 rho_R tau alpha rr k rho_q rho_A rho_ys rho_pies;
-
-psi1 = 1.54;
-psi2 = 0.25;
-psi3 = 0.25;
-rho_R = 0.5;
-alpha = 0.3;
-rr = 2.51;
-k = 0.5;
-tau = 0.5;
-rho_q = 0.4;
-rho_A = 0.2;
-rho_ys = 0.9;
-rho_pies = 0.7;
-
-
-model(linear);
-y = y(+1) - (tau +alpha*(2-alpha)*(1-tau))*(R-pie(+1))-alpha*(tau +alpha*(2-alpha)*(1-tau))*dq(+1) + alpha*(2-alpha)*((1-tau)/tau)*(y_s-y_s(+1))-A(+1);
-pie = exp(-rr/400)*pie(+1)+alpha*exp(-rr/400)*dq(+1)-alpha*dq+(k/(tau+alpha*(2-alpha)*(1-tau)))*y+k*alpha*(2-alpha)*(1-tau)/(tau*(tau+alpha*(2-alpha)*(1-tau)))*y_s;
-pie = de+(1-alpha)*dq+pie_s;
-R = rho_R*R(-1)+(1-rho_R)*(psi1*pie+psi2*(y+alpha*(2-alpha)*((1-tau)/tau)*y_s)+psi3*de)+e_R;
-dq = rho_q*dq(-1)+e_q;
-y_s = rho_ys*y_s(-1)+e_ys;
-pie_s = rho_pies*pie_s(-1)+e_pies;
-A = rho_A*A(-1)+e_A;
-y_obs = y-y(-1)+A;
-pie_obs = 4*pie;
-R_obs = 4*R;
-end;
-
-shocks;
-	var e_R = 1.25^2;
-	var e_q = 2.5^2;
-	var e_A = 1.89;
-	var e_ys = 1.89;
-	var e_pies = 1.89;
-end;
-
-varobs y_obs R_obs pie_obs dq de;
+@#include "ls2003_model.inc"
 
 estimated_params;
 psi1 , gamma_pdf,1.5,0.5;
diff --git a/tests/gsa/ls2003ide.mod b/tests/gsa/ls2003ide.mod
index 134e2240cfae5b6f2d0f0930ee7a3cffe5475e2b..6957cd8d63c56fe4a9be9be1fce70fa1b502f84d 100644
--- a/tests/gsa/ls2003ide.mod
+++ b/tests/gsa/ls2003ide.mod
@@ -1,45 +1,4 @@
-var y y_s R pie dq pie_s de A y_obs pie_obs R_obs;
-varexo e_R e_q e_ys e_pies e_A;
-
-parameters psi1 psi2 psi3 rho_R tau alpha rr k rho_q rho_A rho_ys rho_pies;
-
-psi1 = 1.54;
-psi2 = 0.25;
-psi3 = 0.25;
-rho_R = 0.5;
-alpha = 0.3;
-rr = 2.51;
-k = 0.5;
-tau = 0.5;
-rho_q = 0.4;
-rho_A = 0.2;
-rho_ys = 0.9;
-rho_pies = 0.7;
-
-
-model(linear);
-y = y(+1) - (tau +alpha*(2-alpha)*(1-tau))*(R-pie(+1))-alpha*(tau +alpha*(2-alpha)*(1-tau))*dq(+1) + alpha*(2-alpha)*((1-tau)/tau)*(y_s-y_s(+1))-A(+1);
-pie = exp(-rr/400)*pie(+1)+alpha*exp(-rr/400)*dq(+1)-alpha*dq+(k/(tau+alpha*(2-alpha)*(1-tau)))*y+k*alpha*(2-alpha)*(1-tau)/(tau*(tau+alpha*(2-alpha)*(1-tau)))*y_s;
-pie = de+(1-alpha)*dq+pie_s;
-R = rho_R*R(-1)+(1-rho_R)*(psi1*pie+psi2*(y+alpha*(2-alpha)*((1-tau)/tau)*y_s)+psi3*de)+e_R;
-dq = rho_q*dq(-1)+e_q;
-y_s = rho_ys*y_s(-1)+e_ys;
-pie_s = rho_pies*pie_s(-1)+e_pies;
-A = rho_A*A(-1)+e_A;
-y_obs = y-y(-1)+A;
-pie_obs = 4*pie;
-R_obs = 4*R;
-end;
-
-shocks;
-	var e_R = 1.25^2;
-	var e_q = 2.5^2;
-	var e_A = 1.89;
-	var e_ys = 1.89;
-	var e_pies = 1.89;
-end;
-
-varobs y_obs R_obs pie_obs dq de;
+@#include "ls2003_model.inc"
 
 estimated_params;
 psi1 , gamma_pdf,1.5,0.5;
@@ -65,15 +24,14 @@ end;
 disp('CREATE SCREENING SAMPLE, CHECK FOR STABILITY AND PERFORM A SCREENING FOR IDENTIFICATION ANALYSIS');
 disp('TYPE II ERRORS')
 disp(' ')
-disp('PRESS ENTER TO CONTUNUE');
-pause(5);
+% disp('PRESS ENTER TO CONTUNUE');
+% pause(5);
 
 dynare_sensitivity(identification=1, morris_nliv=6, morris_ntra=50);
 
 disp('CREATE MC SAMPLE, CHECK FOR STABILITY AND PERFORM IDENTIFICATION ANALYSIS');
-disp('WIth analytic derivatives')
+disp('With analytic derivatives')
 disp(' ')
-disp('PRESS ENTER TO CONTUNUE');
-pause(5);
-dynare_sensitivity(identification=1, morris=2);
-
+% disp('PRESS ENTER TO CONTUNUE');
+% pause(5);
+dynare_sensitivity(identification=1, morris=2);
\ No newline at end of file
diff --git a/tests/gsa/ls2003scr.mod b/tests/gsa/ls2003scr.mod
index 611c9753b989727c04a6214ee92a877f17614a2a..3811889efca7405bc461eaa3ae36ca1743404d8c 100644
--- a/tests/gsa/ls2003scr.mod
+++ b/tests/gsa/ls2003scr.mod
@@ -1,45 +1,4 @@
-var y y_s R pie dq pie_s de A y_obs pie_obs R_obs;
-varexo e_R e_q e_ys e_pies e_A;
-
-parameters psi1 psi2 psi3 rho_R tau alpha rr k rho_q rho_A rho_ys rho_pies;
-
-psi1 = 1.54;
-psi2 = 0.25;
-psi3 = 0.25;
-rho_R = 0.5;
-alpha = 0.3;
-rr = 2.51;
-k = 0.5;
-tau = 0.5;
-rho_q = 0.4;
-rho_A = 0.2;
-rho_ys = 0.9;
-rho_pies = 0.7;
-
-
-model(linear);
-y = y(+1) - (tau +alpha*(2-alpha)*(1-tau))*(R-pie(+1))-alpha*(tau +alpha*(2-alpha)*(1-tau))*dq(+1) + alpha*(2-alpha)*((1-tau)/tau)*(y_s-y_s(+1))-A(+1);
-pie = exp(-rr/400)*pie(+1)+alpha*exp(-rr/400)*dq(+1)-alpha*dq+(k/(tau+alpha*(2-alpha)*(1-tau)))*y+k*alpha*(2-alpha)*(1-tau)/(tau*(tau+alpha*(2-alpha)*(1-tau)))*y_s;
-pie = de+(1-alpha)*dq+pie_s;
-R = rho_R*R(-1)+(1-rho_R)*(psi1*pie+psi2*(y+alpha*(2-alpha)*((1-tau)/tau)*y_s)+psi3*de)+e_R;
-dq = rho_q*dq(-1)+e_q;
-y_s = rho_ys*y_s(-1)+e_ys;
-pie_s = rho_pies*pie_s(-1)+e_pies;
-A = rho_A*A(-1)+e_A;
-y_obs = y-y(-1)+A;
-pie_obs = 4*pie;
-R_obs = 4*R;
-end;
-
-shocks;
-	var e_R = 1.25^2;
-	var e_q = 2.5^2;
-	var e_A = 1.89;
-	var e_ys = 1.89;
-	var e_pies = 1.89;
-end;
-
-varobs y_obs R_obs pie_obs dq de;
+@#include "ls2003_model.inc"
 
 estimated_params;
 psi1 , gamma_pdf,1.5,0.5;
@@ -63,8 +22,8 @@ end;
 
   
 disp('CREATE SCREENING SAMPLE, CHECK FOR STABILITY AND PERFORM SENSITIVITY ANALYSIS');
-disp('PRESS ENTER TO CONTUNUE');
-pause;
+% disp('PRESS ENTER TO CONTUNUE');
+% pause;
 
 dynare_sensitivity(morris=1, morris_nliv=6, morris_ntra=20, redform=1, 
   namendo=(:), namexo=(:), namlagendo=(:));