proc iml;
reset log;
reset noname;
/*THREE-STEP AND ITERATIVE STOCHASTIC TEST FOR WEAK SEPARABILITY
				Philippe de Peretti, version 1.0
					peretti@univ-paris1.fr
		STEP 1 : AJUSTMENT FOR SUB-UTILITY SUBJECT,
		STEP 2 : AJUSTMENT FOR OVERALL UTILITY, GIVING R²=R²*,
additional constraints for step 1, sub-budget constraint
		STEP 3 : AJUSTMENT FOR SEPARABILITY SUBJECT TO R=R* 
		(X1,P1) : OTHER GOODS
		(X2,P2) : GOODS WEAKLY SEPARABLE								*/

start WS(X1,X2,P1,P2);
XX1=X1;XX2=X2;X=X1||X2;P=P1||P2;XX=X;
load module=Q;
					/*SUB-ROUTINE 1 : QUADRATIC PROGRAMMING*/
start qp(noms,c,h,g,rel,b,teta);
reset noprint;
st=0.1;
if min(eigval(h))<0 then do;
print 'erreur  la valeur propre miminale est negative';
print 'H n"est pas semidefinie positive';
print 'le programme d"optimisation s"arrete';
stop;
end;
nr=nrow(g);
nc=ncol(g);
                /*mise en forme cannonique*/
rev=(rel='<=');
adj=(-1*rev)+^rev;
g=adj#g;
b=adj#b;
eq=( rel = '=' );
if max(eq)=1 then do;
g=g // -(diag(eq)*g)[loc(eq),];
b=b // -(diag(eq)*b)[loc(eq)];
end;

m=(h || -g`) // (g || j(nrow(g),nrow(g),0));
q=c // -b;

/*resoudre le problème*/
call lcp(rc,w,z,M,q);

/*reporter la solution*/
/*reset noname;
print ({              'LA SOLUTION EST OPTIMALE',
                      'PAS DE SOLUTION POSSIBLE',
                      '',
                      '',
                      '',
                      'LA SOLUTION EST NUMERIQUEMENT INSTABLE',
                      'PAS ASSEZ DE MEMOIRE',
                      'LE NOMBRE D"ITERATIONS EXCEDE'}[rc+1]);*/
reset noname;
teta=z[1:nc];
cste=ncol(G);
objval=(c`*teta)+((teta`*H*teta)/2)+cste;
p=probchi((abs(objval/(st*st))),cste);
seuil=1-p;
/*print 'type de fonction objectif c`*teta+((teta/X)-1)`*H*((teta/X)-1)';
print ,'valeur de la fonction objectif:'objval,
       'valeur du Chi² pour'cste'degres de liberte :'p,
       'p(x<X)='seuil,
        'variable de decision' teta /*[r=noms],
        '**********************************';*/
store seuil objval teta;
finish qp;
					/*SUB-ROUTINE 2 : IID TESTS*/
reset noname;
start IID(X1j);
tau=int(min(nrow(X1j)/2-2,40));
tau=max(tau,4);
			/*Independence test first order: Ljung-Box (1978)*/
rho=j(tau,1,0);
stat=j(tau,1,0);
stat1=j(tau,1,0);
pstat1=j(tau,1,0);
do o=1 to tau;
rho[o]=((X1j[o+1:nrow(X1j)]-(sum(X1j)/nrow(X1j)))`*(X1j[1:nrow(X1j)-o]-(sum(X1j)/nrow(X1j))))/
((X1j-(sum(X1j)/nrow(X1j)))`*(X1j-(sum(X1j)/nrow(X1j))));
stat[o]=rho[o]##2/(nrow(X1j)-o);
stat1[o]=sum(stat[1:o])*nrow(X1j)*(nrow(X1j)+2);
pstat1[o]=1-probchi(stat1[o],o);
end;
ljung=stat1||pstat1;

maljung=ljung[1:4,]//ljung[tau,];
a=j(4,1,'tau=');b=(1:4)`;
rmaljung=concat(a,char(b,4))//concat('tau=',char(tau,4));
Cmaljung={'Q-Stat','P-value','ML-Stat','P-value'};

ljung=ljung[1:4,];

			/*Independence test : second order McLeod-Li (1983)*/
X2j=(X1j-sum(X1j)/nrow(X1j))##2;
rho=j(tau,1,0);
stat=j(tau,1,0);
stat1=j(tau,1,0);
pstat1=j(tau,1,0);
do o=1 to tau;
rho[o]=((X2j[o+1:nrow(X2j)]-(sum(X2j)/nrow(X2j)))`*(X2j[1:nrow(X2j)-o]-(sum(X2j)/nrow(X2j))))/
((X2j-(sum(X2j)/nrow(X2j)))`*(X2j-(sum(X2j)/nrow(X2j))));
stat[o]=rho[o]##2/(nrow(X2j)-o);
stat1[o]=sum(stat[1:o])*nrow(X2j)*(nrow(X2j)+2);
pstat1[o]=1-probchi(stat1[o],o);
end;
McLeod=stat1||pstat1;

maleod=McLeod[1:4,]//McLeod[tau,];
MbT=maljung||maleod;

McLeod=McLeod[1:4,];

					/*TEST IID BASED ON F-STAT*/
/*heterogeneite du premier ordre AR1+trend*/
AR1=X1j[2:nrow(X1j)];
YAR1=X1j[1:nrow(X1j)-1]||(1:nrow(X1j)-1)`||j(nrow(AR1),1,1);

betaAR1=inv(YAR1`*YAR1)*(YAR1`*AR1);
residAR1=AR1-YAR1*betaAR1;

YAR1b=j(nrow(AR1),1,1);
betaAR1b=inv(YAR1b`*YAR1b)*(YAR1b`*AR1);
residAR1b=AR1-YAR1b*betaAR1b;

WaldAR1=(nrow(residAR1b)-ncol(YAR1))*(residAR1b`*residAR1b-residAR1`*residAR1)/(residAR1`*residAR1);
FisherAR1=((residAR1b`*residAR1b-residAR1`*residAR1)/2)/
		((residAR1`*residAR1)/(nrow(residAR1b)-ncol(YAR1)));
PfishAR1=1-probf(fisherAR1,2,nrow(residAR1b)-ncol(YAR1));
PwaldAR1=1-probchi(WaldAR1,2);

if nrow(X1j)>=20 then do;
/*heterogeneite du premier ordre AR2+trend*/
AR2=X1j[3:nrow(X1j)];
YAR2=X1j[1:nrow(X1j)-2]||X1j[2:nrow(X1j)-1]||(1:nrow(X1j)-2)`||j(nrow(AR2),1,1);
betaAR2=inv(YAR2`*YAR2)*(YAR2`*AR2);
residAR2=AR2-YAR2*betaAR2;

YAR2b=j(nrow(AR2),1,1);
betaAR2b=inv(YAR2b`*YAR2b)*(YAR2b`*AR2);
residAR2b=AR2-YAR2b*betaAR2b;

WaldAR2=(nrow(residAR2b)-ncol(YAR2))*(residAR2b`*residAR2b-residAR2`*residAR2)/(residAR2`*residAR2);
FisherAR2=((residAR2b`*residAR2b-residAR2`*residAR2)/3)/
		((residAR2`*residAR2)/(nrow(residAR2b)-ncol(YAR2)));
PfishAR2=1-probf(fisherAR2,3,nrow(residAR2b)-ncol(YAR2));
PwaldAR2=1-probchi(WaldAR2,3);
end;
else do;
WaldAR2=.;
FisherAR2=.;
PfishAR2=.;
PwaldAR2=.;
end;

if nrow(X1j)>=30 then do;
/*heterogeneite du premier ordre AR3+trend*/
AR3=X1j[4:nrow(X1j)];
YAR3=X1j[1:nrow(X1j)-3]||X1j[2:nrow(X1j)-2]||X1j[3:nrow(X1j)-1]||(1:nrow(X1j)-3)`||j(nrow(AR3),1,1);
betaAR3=inv(YAR3`*YAR3)*(YAR3`*AR3);
residAR3=AR3-YAR3*betaAR3;

YAR3b=j(nrow(AR3),1,1);
betaAR3b=inv(YAR3b`*YAR3b)*(YAR3b`*AR3);
residAR3b=AR3-YAR3b*betaAR3b;

WaldAR3=(nrow(residAR3b)-ncol(YAR3))*(residAR3b`*residAR3b-residAR3`*residAR3)/(residAR3`*residAR3);
FisherAR3=((residAR3b`*residAR3b-residAR3`*residAR3)/4)/
		((residAR3`*residAR3)/(nrow(residAR3b)-ncol(YAR3)));
PfishAR3=1-probf(fisherAR3,4,nrow(residAR3b)-ncol(YAR3));
PwaldAR3=1-probchi(WaldAR3,4);
end;
else do;
WaldAR3=.;
FisherAR3=.;
PfishAR3=.;
PwaldAR3=.;
end;

if nrow(X1j)>=40 then do;
/*heterogeneite du premier ordre AR4+trend*/
AR4=X1j[5:nrow(X1j)];
YAR4=X1j[1:nrow(X1j)-4]||X1j[2:nrow(X1j)-3]||X1j[3:nrow(X1j)-2]||X1j[4:nrow(X1j)-1]||(1:nrow(X1j)-4)`||j(nrow(AR4),1,1);
betaAR4=inv(YAR4`*YAR4)*(YAR4`*AR4);
residAR4=AR4-YAR4*betaAR4;

YAR4b=j(nrow(AR4),1,1);
betaAR4b=inv(YAR4b`*YAR4b)*(YAR4b`*AR4);
residAR4b=AR4-YAR4b*betaAR4b;

WaldAR4=(nrow(residAR4b)-ncol(YAR4))*(residAR4b`*residAR4b-residAR4`*residAR4)/(residAR4`*residAR4);
FisherAR4=((residAR4b`*residAR4b-residAR4`*residAR4)/5)/
		((residAR4`*residAR4)/(nrow(residAR4b)-ncol(YAR4)));
PfishAR4=1-probf(fisherAR4,5,nrow(residAR4b)-ncol(YAR4));
PwaldAR4=1-probchi(WaldAR4,5);
end;
else do;
WaldAR4=.;
FisherAR4=.;
PfishAR4=.;
PwaldAR4=.;
end;

rhetero={'F-stat','P-value','Wald','P-value'};
chetero={'1st order, AR1 & trend :','1st order, AR2 & trend :','1st
order, AR3 & trend :','1st order, AR4 & trend :'};
chetero1={'2nd order, AR1 & trend :','2nd order, AR2 & trend :','2nd order, AR3 & trend :','2nd order, AR4 & trend :'};
hetero1=(FisherAR1//pfishAR1//WaldAR1//PWaldAR1)||
(FisherAR2//pfishAR2//WaldAR2//PWaldAR2)||
(FisherAR3//pfishAR3//WaldAR3//PWaldAR3)||
(Fisherar4//pfishAR4//WaldAR4//PWaldAR4);

if nrow(X1j)>=40 then do;
/*séléction de modèle AR4 contre AR3*/
AR4=X1j[5:nrow(X1j)];
YAR4=X1j[1:nrow(X1j)-4]||X1j[2:nrow(X1j)-3]||X1j[3:nrow(X1j)-2]||X1j[4:nrow(X1j)-1]||(1:nrow(X1j)-4)`||j(nrow(AR4),1,1);
betaAR4=inv(YAR4`*YAR4)*(YAR4`*AR4);
residAR4=AR4-YAR4*betaAR4;

YAR4b=X1j[2:nrow(X1j)-3]||X1j[3:nrow(X1j)-2]||X1j[4:nrow(X1j)-1]||(1:nrow(X1j)-4)`||j(nrow(AR4),1,1);
betaAR4b=inv(YAR4b`*YAR4b)*(YAR4b`*AR4);
residAR4b=AR4-YAR4b*betaAR4b;

WaldAR4=(nrow(residAR4b)-ncol(YAR4))*(residAR4b`*residAR4b-residAR4`*residAR4)/(residAR4`*residAR4);
FisherAR4=((residAR4b`*residAR4b-residAR4`*residAR4)/1)/
		((residAR4`*residAR4)/(nrow(residAR4b)-ncol(YAR4)));
PfishAR4=1-probf(fisherAR4,1,nrow(residAR4b)-ncol(YAR4));
PwaldAR4=1-probchi(WaldAR4,1);
end;
else do;
WaldAR4=.;
FisherAR4=.;
PfishAR4=.;
PwaldAR4=.;
end;


if nrow(X1j)>=30 then do;
/*séléction de modèle AR3 contre AR2*/
AR3=X1j[4:nrow(X1j)];
YAR3=X1j[1:nrow(X1j)-3]||X1j[2:nrow(X1j)-2]||X1j[3:nrow(X1j)-1]||(1:nrow(X1j)-3)`||j(nrow(AR3),1,1);
betaAR3=inv(YAR3`*YAR3)*(YAR3`*AR3);
residAR3=AR3-YAR3*betaAR3;

YAR3b=X1j[2:nrow(X1j)-2]||X1j[3:nrow(X1j)-1]||(1:nrow(X1j)-3)`||j(nrow(AR3),1,1);
betaAR3b=inv(YAR3b`*YAR3b)*(YAR3b`*AR3);
residAR3b=AR3-YAR3b*betaAR3b;

WaldAR3=(nrow(residAR3b)-ncol(YAR3))*(residAR3b`*residAR3b-residAR3`*residAR3)/(residAR3`*residAR3);
FisherAR3=((residAR3b`*residAR3b-residAR3`*residAR3)/1)/
		((residAR3`*residAR3)/(nrow(residAR3b)-ncol(YAR3)));
PfishAR3=1-probf(fisherAR3,1,nrow(residAR3b)-ncol(YAR3));
PwaldAR3=1-probchi(WaldAR3,1);
end;
else do;
WaldAR3=.;
FisherAR3=.;
PfishAR3=.;
PwaldAR3=.;
end;

if nrow(X1j)>=20 then do;
/*séléction de modèle AR2 contre AR1*/
AR2=X1j[3:nrow(X1j)];
YAR2=X1j[1:nrow(X1j)-2]||X1j[2:nrow(X1j)-1]||(1:nrow(X1j)-2)`||j(nrow(AR2),1,1);
betaAR2=inv(YAR2`*YAR2)*(YAR2`*AR2);
residAR2=AR2-YAR2*betaAR2;

YAR2b=X1j[2:nrow(X1j)-1]||(1:nrow(X1j)-2)`||j(nrow(AR2),1,1);
betaAR2b=inv(YAR2b`*YAR2b)*(YAR2b`*AR2);
residAR2b=AR2-YAR2b*betaAR2b;

WaldAR2=(nrow(residAR2b)-ncol(YAR2))*(residAR2b`*residAR2b-residAR2`*residAR2)/(residAR2`*residAR2);
FisherAR2=((residAR2b`*residAR2b-residAR2`*residAR2)/(ncol(YAR2)-ncol(YAR2b)))/
		((residAR2`*residAR2)/(nrow(residAR2b)-ncol(YAR2)));
PfishAR2=1-probf(fisherAR2,(ncol(YAR2)-ncol(YAR2b)),nrow(residAR2b)-ncol(YAR2));
PwaldAR2=1-probchi(WaldAR2,(ncol(YAR2)-ncol(YAR2b)));
end;
else do;
WaldAR2=.;
FisherAR2=.;
PfishAR2=.;
PwaldAR2=.;
end;

rshetero={'F-stat','P-value','Wald','P-value'};
cshetero={' Ho:AR4=AR3,Ha:AR4 ',' Ho:AR3=AR2,Ha:AR3, ',' Ho:AR2=AR1,Ha:AR2, '};
shetero1=(Fisherar4//pfishAR4//WaldAR4//PWaldAR4)||
(FisherAR3//pfishAR3//WaldAR3//PWaldAR3)||
(FisherAR2//pfishAR2//WaldAR2//PWaldAR2);

/*heterogeneite du second ordre AR1+trend*/
AR1=X1j[2:nrow(X1j)]##2;
YAR1=X1j[1:nrow(X1j)-1]##2||(1:nrow(X1j)-1)`||j(nrow(AR1),1,1);

betaAR1=inv(YAR1`*YAR1)*(YAR1`*AR1);
residAR1=AR1-YAR1*betaAR1;

YAR1b=j(nrow(AR1),1,1);
betaAR1b=inv(YAR1b`*YAR1b)*(YAR1b`*AR1);
residAR1b=AR1-YAR1b*betaAR1b;

WaldAR1=(nrow(residAR1b)-ncol(YAR1))*(residAR1b`*residAR1b-residAR1`*residAR1)/(residAR1`*residAR1);
FisherAR1=((residAR1b`*residAR1b-residAR1`*residAR1)/2)/
		((residAR1`*residAR1)/(nrow(residAR1b)-ncol(YAR1)));
PfishAR1=1-probf(fisherAR1,2,nrow(residAR1b)-ncol(YAR1));
PwaldAR1=1-probchi(WaldAR1,2);

if nrow(X1j)>=20 then do;
/*heterogeneite du second ordre AR2+trend*/
AR2=X1j[3:nrow(X1j)]##2;
YAR2=X1j[1:nrow(X1j)-2]##2||X1j[2:nrow(X1j)-1]##2||X1j[1:nrow(X1j)-2]#X1j[2:nrow(X1j)-1]||
(1:nrow(X1j)-2)`||j(nrow(AR2),1,1);
betaAR2=inv(YAR2`*YAR2)*(YAR2`*AR2);
residAR2=AR2-YAR2*betaAR2;

YAR2b=j(nrow(AR2),1,1);
betaAR2b=inv(YAR2b`*YAR2b)*(YAR2b`*AR2);
residAR2b=AR2-YAR2b*betaAR2b;

WaldAR2=(nrow(residAR2b)-ncol(YAR2))*(residAR2b`*residAR2b-residAR2`*residAR2)/(residAR2`*residAR2);
FisherAR2=((residAR2b`*residAR2b-residAR2`*residAR2)/4)/
		((residAR2`*residAR2)/(nrow(residAR2b)-ncol(YAR2)));
PfishAR2=1-probf(fisherAR2,4,nrow(residAR2b)-ncol(YAR2));
PwaldAR2=1-probchi(WaldAR2,4);
end;
else do;
WaldAR2=.;
FisherAR2=.;
PfishAR2=.;
PwaldAR2=.;
end;


if nrow(X1j)>=30 then do;
/*heterogeneite du second ordre AR3+trend*/
AR3=X1j[4:nrow(X1j)]##2;
YAR3=X1j[1:nrow(X1j)-3]##2||X1j[2:nrow(X1j)-2]##2||X1j[3:nrow(X1j)-1]##2||
X1j[1:nrow(X1j)-3]#X1j[2:nrow(X1j)-2]||X1j[1:nrow(X1j)-3]#X1j[3:nrow(X1j)-1]||
X1j[2:nrow(X1j)-2]#X1j[3:nrow(X1j)-1]
||(1:nrow(X1j)-3)`||j(nrow(AR3),1,1);
betaAR3=inv(YAR3`*YAR3)*(YAR3`*AR3);
residAR3=AR3-YAR3*betaAR3;

YAR3b=j(nrow(AR3),1,1);
betaAR3b=inv(YAR3b`*YAR3b)*(YAR3b`*AR3);
residAR3b=AR3-YAR3b*betaAR3b;

WaldAR3=(nrow(residAR3b)-ncol(YAR3))*(residAR3b`*residAR3b-residAR3`*residAR3)/(residAR3`*residAR3);
FisherAR3=((residAR3b`*residAR3b-residAR3`*residAR3)/7)/
		((residAR3`*residAR3)/(nrow(residAR3b)-ncol(YAR3)));
PfishAR3=1-probf(fisherAR3,7,nrow(residAR3b)-ncol(YAR3));
PwaldAR3=1-probchi(WaldAR3,7);
end;
else do;
WaldAR3=.;
FisherAR3=.;
PfishAR3=.;
PwaldAR3=.;
end;


if nrow(X1j)>=40 then do;
/*heterogeneite du second ordre AR4+trend*/
AR4=X1j[5:nrow(X1j)]##2;
YAR4=X1j[1:nrow(X1j)-4]##2||X1j[2:nrow(X1j)-3]##2||X1j[3:nrow(X1j)-2]##2||X1j[4:nrow(X1j)-1]##2||
X1j[4:nrow(X1j)-1]#X1j[3:nrow(X1j)-2]||X1j[4:nrow(X1j)-1]#X1j[2:nrow(X1j)-3]||X1j[4:nrow(X1j)-1]#X1j[1:nrow(X1j)-4]||
X1j[3:nrow(X1j)-2]#X1j[2:nrow(X1j)-3]||X1j[3:nrow(X1j)-2]#X1j[1:nrow(X1j)-4]||
X1j[1:nrow(X1j)-4]#X1j[2:nrow(X1j)-3]||
(1:nrow(X1j)-4)`||j(nrow(AR4),1,1);

betaAR4=inv(YAR4`*YAR4)*(YAR4`*AR4);
residAR4=AR4-YAR4*betaAR4;

YAR4b=j(nrow(AR4),1,1);
betaAR4b=inv(YAR4b`*YAR4b)*(YAR4b`*AR4);
residAR4b=AR4-YAR4b*betaAR4b;

WaldAR4=(nrow(residAR4b)-ncol(YAR4))*(residAR4b`*residAR4b-residAR4`*residAR4)/(residAR4`*residAR4);
FisherAR4=((residAR4b`*residAR4b-residAR4`*residAR4)/11)/
		((residAR4`*residAR4)/(nrow(residAR4b)-ncol(YAR4)));
PfishAR4=1-probf(fisherAR4,11,nrow(residAR4b)-ncol(YAR4));
PwaldAR4=1-probchi(WaldAR4,11);
end;
else do;
WaldAR4=.;
FisherAR4=.;
PfishAR4=.;
PwaldAR4=.;
end;



hetero2=(FisherAR1//pfishAR1//WaldAR1//PWaldAR1)||
(FisherAR2//pfishAR2//WaldAR2//PWaldAR2)||
(FisherAR3//pfishAR3//WaldAR3//PWaldAR3)||
(Fisherar4//pfishAR4//WaldAR4//PWaldAR4);

if nrow(X1j)>=40 then do;
/*selection des modèles AR4 contre AR3*/
AR4=X1j[5:nrow(X1j)]##2;
YAR4=X1j[1:nrow(X1j)-4]##2||X1j[2:nrow(X1j)-3]##2||X1j[3:nrow(X1j)-2]##2||X1j[4:nrow(X1j)-1]##2||
X1j[4:nrow(X1j)-1]#X1j[3:nrow(X1j)-2]||X1j[4:nrow(X1j)-1]#X1j[2:nrow(X1j)-3]||X1j[4:nrow(X1j)-1]#X1j[1:nrow(X1j)-4]||
X1j[3:nrow(X1j)-2]#X1j[2:nrow(X1j)-3]||X1j[3:nrow(X1j)-2]#X1j[1:nrow(X1j)-4]||
X1j[1:nrow(X1j)-4]#X1j[2:nrow(X1j)-3]||
(1:nrow(X1j)-4)`||j(nrow(AR4),1,1);

betaAR4=inv(YAR4`*YAR4)*(YAR4`*AR4);
residAR4=AR4-YAR4*betaAR4;

YAR4b=X1j[2:nrow(X1j)-3]##2||X1j[3:nrow(X1j)-2]##2||X1j[4:nrow(X1j)-1]##2||
X1j[4:nrow(X1j)-1]#X1j[3:nrow(X1j)-2]||X1j[4:nrow(X1j)-1]#X1j[2:nrow(X1j)-3]||
X1j[3:nrow(X1j)-2]#X1j[2:nrow(X1j)-3]||
(1:nrow(X1j)-4)`||j(nrow(AR4),1,1);

betaAR4b=inv(YAR4b`*YAR4b)*(YAR4b`*AR4);
residAR4b=AR4-YAR4b*betaAR4b;

WaldAR4=(nrow(residAR4b)-ncol(YAR4))*(residAR4b`*residAR4b-residAR4`*residAR4)/(residAR4`*residAR4);
FisherAR4=((residAR4b`*residAR4b-residAR4`*residAR4)/(ncol(YAR4)-ncol(YAR4B))/
		((residAR4`*residAR4)/(nrow(residAR4b)-ncol(YAR4))));
PfishAR4=1-probf(fisherAR4,(ncol(YAR4)-ncol(YAR4B)),nrow(residAR4b)-ncol(YAR4));
PwaldAR4=1-probchi(WaldAR4,(ncol(YAR4)-ncol(YAR4B)));
end;
else do;
WaldAR4=.;
FisherAR4=.;
PfishAR4=.;
PwaldAR4=.;
end;


if nrow(X1j)>=30 then do;
/*selection des modèles AR3 contre AR2*/
AR3=X1j[4:nrow(X1j)]##2;
YAR3=X1j[1:nrow(X1j)-3]##2||X1j[2:nrow(X1j)-2]##2||X1j[3:nrow(X1j)-1]##2||
X1j[1:nrow(X1j)-3]#X1j[2:nrow(X1j)-2]||X1j[1:nrow(X1j)-3]#X1j[3:nrow(X1j)-1]||
X1j[2:nrow(X1j)-2]#X1j[3:nrow(X1j)-1]
||(1:nrow(X1j)-3)`||j(nrow(AR3),1,1);
betaAR3=inv(YAR3`*YAR3)*(YAR3`*AR3);
residAR3=AR3-YAR3*betaAR3;

YAR3b=X1j[2:nrow(X1j)-2]##2||X1j[3:nrow(X1j)-1]##2||
X1j[2:nrow(X1j)-2]#X1j[3:nrow(X1j)-1]
||(1:nrow(X1j)-3)`||j(nrow(AR3),1,1);
betaAR3b=inv(YAR3b`*YAR3b)*(YAR3b`*AR3);
residAR3b=AR3-YAR3b*betaAR3b;

WaldAR3=(nrow(residAR3b)-ncol(YAR3))*(residAR3b`*residAR3b-residAR3`*residAR3)/(residAR3`*residAR3);
FisherAR3=((residAR3b`*residAR3b-residAR3`*residAR3)/(ncol(YAR3)-ncol(YAR3b)))/
		((residAR3`*residAR3)/(nrow(residAR3b)-ncol(YAR3)));
PfishAR3=1-probf(fisherAR3,ncol(YAR3)-ncol(YAR3b),nrow(residAR3b)-ncol(YAR3));
PwaldAR3=1-probchi(WaldAR3,ncol(YAR3)-ncol(YAR3b));
end;
else do;
WaldAR3=.;
FisherAR3=.;
PfishAR3=.;
PwaldAR3=.;
end;

if nrow(X1j)>=20 then do;
/*selection des modèles AR2 contre AR1*/
AR2=X1j[3:nrow(X1j)]##2;
YAR2=X1j[1:nrow(X1j)-2]##2||X1j[2:nrow(X1j)-1]##2||X1j[1:nrow(X1j)-2]#X1j[2:nrow(X1j)-1]||
(1:nrow(X1j)-2)`||j(nrow(AR2),1,1);
betaAR2=inv(YAR2`*YAR2)*(YAR2`*AR2);
residAR2=AR2-YAR2*betaAR2;

YAR2b=X1j[2:nrow(X1j)-1]##2||(1:nrow(X1j)-2)`||j(nrow(AR2),1,1);
betaAR2b=inv(YAR2b`*YAR2b)*(YAR2b`*AR2);
residAR2b=AR2-YAR2b*betaAR2b;

WaldAR2=(nrow(residAR2b)-ncol(YAR2))*(residAR2b`*residAR2b-residAR2`*residAR2)/(residAR2`*residAR2);
FisherAR2=((residAR2b`*residAR2b-residAR2`*residAR2)/(ncol(YAR2)-ncol(YAR2b)))/
		((residAR2`*residAR2)/(nrow(residAR2b)-ncol(YAR2)));
PfishAR2=1-probf(fisherAR2,(ncol(YAR2)-ncol(YAR2b)),nrow(residAR2b)-ncol(YAR2));
PwaldAR2=1-probchi(WaldAR2,(ncol(YAR2)-ncol(YAR2b)));
end;
else do;
WaldAR2=.;
FisherAR2=.;
PfishAR2=.;
PwaldAR2=.;
end;
shetero2=(Fisherar4//pfishAR4//WaldAR4//PWaldAR4)||
(FisherAR3//pfishAR3//WaldAR3//PWaldAR3)||
(FisherAR2//pfishAR2//WaldAR2//PWaldAR2);
hetero1=hetero1`;hetero2=hetero2`;
MbT=maljung||maleod;

/*selection automatique des retards pour les regressions auxiliaires*/
tau1=4;tau2=4;
do i=1 to nrow(shetero1);
	do j=1 to ncol(shetero1);
	if shetero1[i,j]=. then shetero1[i,j]=1;
	else shetero1[i,j]=shetero1[i,j];
	if shetero2[i,j]=. then shetero2[i,j]=1;
	else shetero2[i,j]=shetero2[i,j];
	end;
end;
do i=1 to ncol(shetero1);
	if max(shetero1[2,i],shetero1[4,i])<0.05 then tau1=tau1;
	else tau1=tau1-1;
	if max(shetero2[2,i],shetero2[4,i])<0.05 then tau2=tau2;
	else tau2=tau2-1;
end;
Fo=hetero1[tau1,]`;rFo=chetero[tau1,];
So=hetero2[tau2,]`;rSo=chetero1[tau2,];
rFo1=rFo||rSo;
Fo1=(Fo||So)`;
tau1=char(tau1,1);
tau2=char(tau2,1);

print 	'STATISTICAL INFERENCE : RESULTS OF IID TESTS',
		'___________________________________________',
		'A-Moment-based tests';
print 	'First (Ljung-Box Q-Stat) and second (McLeod & Li ML-Stat) order',
		MbT[rowname=Rmaljung colname=Cmaljung],
		'___________________________________________';
print 	'B-Auxiliary regressions',
'First order at lag' tau1 ', second order at lag' tau2;
print Fo1[rowname=rFo1 colname=rhetero],
'_________________________________________________________________';
store MbT Fo1;
finish IID;
					/*SUB-ROUTINE 3 : Q-INDEX*/

/*programme de calcul des indices d'utilité et de prix*/
start Q(P,X);
W=j(nrow(P),ncol(P),0);
WB=j(nrow(P),ncol(P),0);
WB1=j(nrow(P),ncol(P),0);
Rev=j(nrow(P),1,0);
Qind=j(nrow(P),1,0);
utilite=j(nrow(P),1,100);
do i=1 to nrow(P);/*calcul du revenu*/
Rev[i]=X[i,]*P[i,]`;
end;
/*calcul des parts budgétaires*/
do i=1 to nrow(P);
do j=1 to ncol(P);
W[i,j]=(X[i,j]*P[i,j])/Rev[i];
end;
end;
do i=2 to nrow(P);/*calcul du revenu*/
do j=1 to ncol(P);
WB[i,j]=(W[i,j]+W[i-1,j])/2;
end;
end;

do i=2 to nrow(P);/*calcul de l'indice*/
do j=1 to ncol(P);
WB1[i,j]=WB[i,j]*(log(X[i,j])-log(X[i-1,j]));
end;
end;

do i=1 to nrow(P);
Qind[i]=sum(WB1[i,]);
end;

do i=2 to nrow(P);
utilite[i]=exp(Qind[i]+log(utilite[i-1]));
end;
R1=(Rev/Rev[1])*100;
Pind=(R1/utilite);
store Pind utilite Qind;
finish q;

					/*STEP 1 : SUB-UTILITY*/

	/*GARP for sub-utility*/
free sol;
SDRP=j(nrow(X2),nrow(X2),0);
DRP=j(nrow(X2),nrow(X2),0);
do i=1 to nrow(X2);
	do j=1 to nrow(X2);
	if (P2[i,]*X2[i,]`) > (P2[i,]*X2[j,]`) then SDRP[i,j]=1;else SDRP[i,j]=0;
	if (P2[i,]*X2[i,]`) >=(P2[i,]*X2[j,]`) then  DRP[i,j]=1;else  DRP[i,j]=0;
	end;
end;
GARP=DRP;
do k=1 to nrow(DRP);
	do i=1 to nrow(DRP);
		do j=1 to nrow(DRP);
		if GARP[i,k]=0 | GARP[k,j]=0 then t=1; else GARP[i,j]=1;
		end;
	end;
end;
	/*# of violations*/
free couple;nvio=0;
do i=1 to nrow(X2);
	do j=1 to nrow(X2);
	if GARP[i,j]=1 & SDRP[j,i]=1 then do;
		nvio=nvio+1;
		couple=couple//(i||j);
		end;
	else t=1;
	end;
end;
nvioS=nvio;
if nvio^=0 then do;
S=1;
/*print couple;*/
end;
else S=0;
do while (nvio^=0);
	/*searching for the maxima*/
	free couple1;
		do i=1 to nrow(couple);
		couple1=couple1//(couple[i,]||sum(GARP[couple[i,1],]));
		end;
	free couple1a;
	m=max(couple1[,3]);
	do i=1 to nrow(couple1);
	if couple1[i,3]=m then couple1a=couple1a//couple1[i,];
	else t=1;
	end;

		/*simplification of couple1*/
	free couple2;
	couple2=couple1a[1,];
	do i=2 to nrow(couple1a);
	if couple1a[i,1]^=couple1a[i-1,1] then couple2=couple2//couple1a[i,];
	else t=1;
	end;
	free mat;
	do i=1 to nrow(couple2);
		c=(-2/XX2[couple2[i,1],])`;
		H=2/(XX2[couple2[i,1],]##2);H=diag(H);
		G=I(ncol(XX2))//P2[couple2[i,1],];
		REL=j(ncol(XX2),1,'>')//'=';
		B=j(ncol(X2),1,0.01)//(P2[couple2[i,1],]*X2[couple2[i,1],]`);
			do j=1 to nrow(GARP);
			if GARP[couple2[i,1],j]=1 & (couple2[i,1]^=j) then do;
				G=G//P2[j,];
				rel=rel//'>=';
				B=B//P2[j,]*X2[j,]`+0.000001;
				end;
			else t=1;
			end;
				/*additional constraints for R=R* */
			/*do j=1 to nrow(GARP);
			if GARPB[couple2[i,1],j]=1 & (couple2[i,1]^=j) then do;
				G=G//P2[j,];
				rel=rel//'>=';
				B=B//(P1[j,]*X1[j,]`+P2[j,]*X2[j,]`-P1[j,]*X1[couple2[i,1],]`+0.000001);
				end;
			else t=1;
			end;*/
		run qp(noms,C,H,G,rel,B,teta);
        load objval teta;
		mat=mat//(objval||couple2[i,]||teta`);
	end;
	mn=min(mat[,1]);
	do i=1 to nrow(mat);
	if mat[i,1]=mn then do;
	X2[mat[i,2],]=mat[i,5:ncol(mat)];
	sol=sol//mat[i,];
	end;
	else t=1;
	end;
	/*GARP with adjusted bundles*/
SDRP=j(nrow(X2),nrow(X2),0);
DRP=j(nrow(X2),nrow(X2),0);
do i=1 to nrow(X2);
	do j=1 to nrow(X2);
	if (P2[i,]*X2[i,]`) > (P2[i,]*X2[j,]`) then SDRP[i,j]=1;else SDRP[i,j]=0;
	if (P2[i,]*X2[i,]`) >=(P2[i,]*X2[j,]`) then  DRP[i,j]=1;else  DRP[i,j]=0;
	end;
end;
GARP=DRP;
do k=1 to nrow(DRP);
	do i=1 to nrow(DRP);
		do j=1 to nrow(DRP);
		if GARP[i,k]=0 | GARP[k,j]=0 then t=1; else GARP[i,j]=1;
		end;
	end;
end;
	/*# of violations*/
free couple;nvio=0;
do i=1 to nrow(X2);
	do j=1 to nrow(X2);
	if GARP[i,j]=1 & SDRP[j,i]=1 then do;
		nvio=nvio+1;
		couple=couple//(i||j);
		end;
	else t=1;
	end;
end;
end;

free adjS outS;
if nvioS^=0 then do;
outS=sol;
do i=1 to nrow(sol);
adjS=adjS//(X2[sol[i,2],]/XX2[sol[i,2],]-1);
end;
store adjS outS;
end;
GARPA=GARP;
X=X1||X2;
P=P1||P2;

					/*STEP 2 : OVERALL UTILITY*/
free sol;
/*1-a : GARP for the overall utility*/
SDRP=j(nrow(X),nrow(X),0);
DRP=j(nrow(X),nrow(X),0);
do i=1 to nrow(X);
	do j=1 to nrow(X);
	if (P[i,]*X[i,]`) > (P[i,]*X[j,]`) then SDRP[i,j]=1;else SDRP[i,j]=0;
	if (P[i,]*X[i,]`) >=(P[i,]*X[j,]`) then  DRP[i,j]=1;else  DRP[i,j]=0;
	end;
end;
GARP=DRP;
do k=1 to nrow(DRP);
	do i=1 to nrow(DRP);
		do j=1 to nrow(DRP);
		if GARP[i,k]=0 | GARP[k,j]=0 then t=1; else GARP[i,j]=1;
		end;
	end;
end;
	/*# of violations*/
free couple;nvio=0;
do i=1 to nrow(X);
	do j=1 to nrow(X);
	if GARP[i,j]=1 & SDRP[j,i]=1 then do;
		nvio=nvio+1;
		couple=couple//(i||j);
		end;
	else t=1;
	end;
end;
nvioU=nvio;
/*if nvio^=0 then do;
print couple;
end;
else t=1;*/
do while (nvio^=0);
	/*searching for the maxima*/
	free couple1;
		do i=1 to nrow(couple);
		couple1=couple1//(couple[i,]||sum(GARP[couple[i,1],]));
		end;
	free couple1a;
	m=max(couple1[,3]);
	do i=1 to nrow(couple1);
	if couple1[i,3]=m then couple1a=couple1a//couple1[i,];
	else t=1;
	end;

	/*simplification of couple1*/
	free couple2;
	couple2=couple1a[1,];
	do i=2 to nrow(couple1a);
	if couple1a[i,1]^=couple1a[i-1,1] then couple2=couple2//couple1a[i,];
	else t=1;
	end;
	free mat;
	do i=1 to nrow(couple2);
		c=(-2/XX[couple2[i,1],])`;
		H=2/(XX[couple2[i,1],]##2);H=diag(H);
		G=I(ncol(XX))//P[couple2[i,1],];
		REL=j(ncol(X),1,'>')//'=';
		B=j(ncol(X),1,0.01)//(P[couple2[i,1],]*X[couple2[i,1],]`);
			do j=1 to nrow(GARP);
			if GARP[couple2[i,1],j]=1 & (couple2[i,1]^=j) then do;
				G=G//P[j,];
				rel=rel//'>=';
				B=B//P[j,]*X[j,]`+0.000001;
				end;
			else t=1;
			end;
		/*additional constraint for sub-budget weakly separable*/
*		G=G//(j(1,ncol(P1),0)||P2[couple2[i,1],]);
*		B=B//(P2[couple2[i,1],]*X2[couple2[i,1],]`);
*		rel=rel//'=';
		/*additional constraint for sub-utility*/
		do j=1 to nrow(GARPA);
			if GARPA[couple2[i,1],j]=1 & (couple2[i,1]^=j) then do;
				G=G//(j(1,ncol(P1),0)||P2[j,]);
				rel=rel//'>=';
				B=B//P2[j,]*X2[j,]`+0.000001;
				end;
			else t=1;
			end;
		run qp(noms,C,H,G,rel,B,teta);
        load objval teta;
		mat=mat//(objval||couple2[i,]||teta`);
	end;
	mn=min(mat[,1]);
	do i=1 to nrow(mat);
	if mat[i,1]=mn then do;
	X[mat[i,2],]=mat[i,5:ncol(mat)];
	sol=sol//mat[i,];
	end;
	else t=1;
	end;
	/*GARP with adjusted bundles*/
SDRP=j(nrow(X),nrow(X),0);
DRP=j(nrow(X),nrow(X),0);
do i=1 to nrow(X);
	do j=1 to nrow(X);
	if (P[i,]*X[i,]`) > (P[i,]*X[j,]`) then SDRP[i,j]=1;else SDRP[i,j]=0;
	if (P[i,]*X[i,]`) >=(P[i,]*X[j,]`) then  DRP[i,j]=1;else  DRP[i,j]=0;
	end;
end;
GARP=DRP;
do k=1 to nrow(DRP);
	do i=1 to nrow(DRP);
		do j=1 to nrow(DRP);
		if GARP[i,k]=0 | GARP[k,j]=0 then t=1; else GARP[i,j]=1;
		end;
	end;
end;
	/*# of violations*/
free couple;nvio=0;
do i=1 to nrow(X);
	do j=1 to nrow(X);
	if GARP[i,j]=1 & SDRP[j,i]=1 then do;
		nvio=nvio+1;
		couple=couple//(i||j);
		end;
	else t=1;
	end;
end;
end;
GARPB=GARP;
free adjU outU;
if nvioU^=0 then do;
outU=sol;
do i=1 to nrow(outU);
adjU=adjU//(X[outU[i,2],]/XX[outU[i,2],]-1);
X1=X[,1:ncol(XX1)];X2=X[,ncol(X1)+1:ncol(X)];
end;

						/*normalisation*/
store outU adjU;
end;
					/*computing utility indices for step-three*/
run Q(P2,X2);
load Pind utilite;
U=utilite;Lambda=Pind;

					/*STEP 3 : SEPARABILITY*/

	/*GARP*/
free sol;
SDRP=j(nrow(X1),nrow(X1),0);
DRP=j(nrow(X1),nrow(X1),0);
do i=1 to nrow(X1);
	do j=1 to nrow(X1);
	if (P1[i,]*X1[i,]`+Lambda[i]*U[i]) > (P1[i,]*X1[j,]`+Lambda[i]*U[j]) then
SDRP[i,j]=1;else SDRP[i,j]=0;
	if (P1[i,]*X1[i,]`+Lambda[i]*U[i]) >=(P1[i,]*X1[j,]`+Lambda[i]*U[j]) then 
DRP[i,j]=1;else  DRP[i,j]=0;
	end;
end;
GARP=DRP;
do k=1 to nrow(DRP);
	do i=1 to nrow(DRP);
		do j=1 to nrow(DRP);
		if GARP[i,k]=0 | GARP[k,j]=0 then t=1; else GARP[i,j]=1;
		end;
	end;
end;
/*# of violations*/
free couple;nvio=0;
do i=1 to nrow(X1);
	do j=1 to nrow(X1);
	if GARP[i,j]=1 & SDRP[j,i]=1 then do;
		nvio=nvio+1;
		couple=couple//(i||j);
		end;
	else t=1;
	end;
end;
nvioSEP=nvio;
if nvio^=0 then do;
SEP=1;
/*print couple;*/
end;
else SEP=0;
do while (nvio^=0);
	/*searching for the maxima*/
	free couple1;
		do i=1 to nrow(couple);
		couple1=couple1//(couple[i,]||sum(GARP[couple[i,1],]));
		end;
	free couple1a;
	m=max(couple1[,3]);
	do i=1 to nrow(couple1);
	if couple1[i,3]=m then couple1a=couple1a//couple1[i,];
	else t=1;
	end;

	/*simplification of couple1*/
	free couple2;
	couple2=couple1a[1,];
	do i=2 to nrow(couple1a);
	if couple1a[i,1]^=couple1a[i-1,1] then couple2=couple2//couple1a[i,];
	else t=1;
	end;
	free mat;
	do i=1 to nrow(couple2);
		c=(-2/XX1[couple2[i,1],])`;
		H=2/(XX1[couple2[i,1],]##2);H=diag(H);
		G=I(ncol(XX1))//P1[couple2[i,1],];
		REL=j(ncol(XX1),1,'>')//'=';
		B=j(ncol(XX1),1,0.01)//(P1[couple2[i,1],]*X1[couple2[i,1],]`);
			do j=1 to nrow(GARP);
			if GARP[couple2[i,1],j]=1 & (couple2[i,1]^=j) then do;
				G=G//P1[j,];
				rel=rel//'>=';
				B=B//(P1[j,]*X1[j,]`+Lambda[j]*U[j]-Lambda[j]*U[couple2[i,1]]+0.000001);
				end;
			else t=1;
			end;
			/*additional constraints for R=R* */
			do j=1 to nrow(GARP);
			if GARPB[couple2[i,1],j]=1 & (couple2[i,1]^=j) then do;
				G=G//P1[j,];
				rel=rel//'>=';
				B=B//(P1[j,]*X1[j,]`+P2[j,]*X2[j,]`-P2[j,]*X2[couple2[i,1],]`+0.000001);
				end;
			else t=1;
			end;
		run qp(noms,C,H,G,rel,B,teta);
        load objval teta;
		mat=mat//(objval||couple2[i,]||teta`);
	end;
	mn=min(mat[,1]);
	do i=1 to nrow(mat);
	if mat[i,1]=mn then do;
	X1[mat[i,2],]=mat[i,5:ncol(mat)];
	sol=sol//mat[i,];
	end;
	else t=1;
	end;
	/*GARP for adjusted bundles*/
SDRP=j(nrow(X1),nrow(X1),0);
DRP=j(nrow(X1),nrow(X1),0);
do i=1 to nrow(X1);
	do j=1 to nrow(X1);
	if (P1[i,]*X1[i,]`+Lambda[i]*U[i]) > (P1[i,]*X1[j,]`+Lambda[i]*U[j]) then
SDRP[i,j]=1;else SDRP[i,j]=0;
	if (P1[i,]*X1[i,]`+Lambda[i]*U[i]) >=(P1[i,]*X1[j,]`+Lambda[i]*U[j]) then 
DRP[i,j]=1;else  DRP[i,j]=0;
	end;
end;
GARP=DRP;
do k=1 to nrow(DRP);
	do i=1 to nrow(DRP);
		do j=1 to nrow(DRP);
		if GARP[i,k]=0 | GARP[k,j]=0 then t=1; else GARP[i,j]=1;
		end;
	end;
end;
	/*# of violations*/
free couple;nvio=0;
do i=1 to nrow(X1);
	do j=1 to nrow(X1);
	if GARP[i,j]=1 & SDRP[j,i]=1 then do;
		nvio=nvio+1;
		couple=couple//(i||j);
		end;
	else t=1;
	end;
end;
end;

free adjSEP outSEP;
if SEP^=0 then do;
outSEP=sol;
do i=1 to nrow(sol);
adjSEP=adjSEP//(X1[sol[i,2],]/XX1[sol[i,2],]-1);
end;

store adjSEP outSEP;
end;
	/*Printings*/
RES1=nvioU||nvioS||nvioSEP;
rRES1={'# of violations'};
cRES1={'Overall utility','Sub-utility','Separability'};
print res1[rowname=rRES1 colname=cRES1];

					/*STATISTICAL ANALYSIS*/
	/*Adjustment for X2*/
if ((nvioU+nvioS)^=0) then do;
	adjX2=(X2/XX2-1);
	free S1_X2 S2_X2;
		do i=1 to ncol(adjX2);
		S1_X2=S1_X2//adjX2[,i];
		end;
		do i=1 to nrow(S1_X2);
		if S1_X2[i]^=0 then S2_X2=S2_X2//S1_X2[i];
		else t=1;
		end;
end;
else do;
S1_X2=j(nrow(X2)*ncol(X2),1,0);
S2_X2=S1_X2;
end;

	/*Adjustment for X1*/
if ((nvioU+nvioSEP)^=0) then do;
	adjX1=(X1/XX1-1);
	free S1_X1 S2_X1;
		do i=1 to ncol(adjX1);
		S1_X1=S1_X1//adjX1[,i];
		end;
		do i=1 to nrow(S1_X1);
		if S1_X1[i]^=0 then S2_X1=S2_X1//S1_X1[i];
		else t=1;
		end;
end;
else do;
S1_X1=j(nrow(X1)*ncol(X1),1,0);
S2_X1=S1_X1;
end;

if nvioU^=0 then outU1=(1:nrow(outU))`||outU[,2]||outU[,1];
else outU1=(.||.);
if nvioS^=0 then outS1=(1:nrow(outS))`||outS[,2]||outS[,1];
else outS1=(.||.);
if nvioSEP^=0 then outSEP1=(1:nrow(outSEP))`||outSEP[,2]||outSEP[,1];
else outSEP1=(.||.);
c01={'Iteration','X*i','Adjustment'};
print 'Overall utility',outU1[colname=c01],
'Sub-utility',outS1[colname=c01],
'Separability',outSEP1[colname=c01];

if ((nvioU+nvioSEP)^=0) then do;
	print 	'STATISTICAL INFERENCE FOR X1';
	run iid(S2_X1);
	load MbT Fo1;
	MbTA=MbT;
	Fo1A=Fo1;
	ADX1=1;
	store MbTA Fo1A;
	end;
	else ADX1=0;
	store ADX1;
if ((nvioU+nvioS)^=0) then do;
	print 	'STATISTICAL INFERENCE FOR X2';
	run iid(S2_X2);
	load MbT Fo1;
	MbTB=MbT;
	Fo1B=Fo1;
	ADX2=1;
	store MbTB Fo1B;
	end;
	else ADX2=0;
	store ADX2;

S1_X1=S1_X1||((S1_X1-sum(S1_X1)/nrow(S1_X1))##2)||(1:nrow(S1_X1))`;
S2_X1=S2_X1||((S2_X1-sum(S2_X1)/nrow(S2_X1))##2)||(1:nrow(S2_X1))`;
S1_X2=S1_X2||((S1_X2-sum(S1_X2)/nrow(S1_X2))##2)||(1:nrow(S1_X2))`;
S2_X2=S2_X2||((S2_X2-sum(S2_X2)/nrow(S2_X2))##2)||(1:nrow(S2_X2))`;

co={'serie','serieQ','date'};

create adjS1X1 from S1_X1[colname=co];
append from S1_X1;
close adjS1X1;

create adjS2X1 from S2_X1[colname=co];
append from S2_X1;
close adjS2X1;

create adjS1X2 from S1_X2[colname=co];
append from S1_X2;
close adjS1X2;

create adjS2X2 from S2_X2[colname=co];
append from S2_X2;
close adjS2X2;

store S2_X1 S2_X2;
finish WS;
store module=(QP WS IID Q);