This code produces the table of results in “On the Non-Identification of Counterfactuals in Dynamic Games” by Kalouptsidi, Scott, and Souza-Rodrigues
please report any errors to ptscott@gmail.com

Model

(* parameters *)
beta=.9; (* discount factor *)
r=.1; (* scrap value *)
c=-.2; (* entry cost *)
piparam1=1.2;(* monopoly profits *)
piparam2=-1.2;(* duopoly profits *)
trans[a_,s_]:=a; (* transition of states *)
var=1;
basetheta={r,c,piparam1,piparam2};

NumericListQ[list_]:=And@@NumericQ/@list;
NumericArrayQ[list_]:=And@@NumericQ/@Flatten[list];

(* convert bivariate state vector to integer *)
s2x[s_]:=Switch[s,
{2,2},1,
{2,1},2,
{1,2},3,
{1,1},4
];

(* convert index of state vector to bivariate *)
x2s[x_]:=Switch[x,
1,{2,2},
2,{2,1},
3,{1,2},
4,{1,1}
];


(* payoff function *)
pii[aown_,aother_,sown_,sother_,theta_,norm_]:=Module[{pay,r,c,pi1,pi2},
r=theta[[1]];
c=theta[[2]];
pi1=theta[[3]];
pi2=theta[[4]];

If[norm==1,
(* PSD version *)
If[aown==2&&sown==2,pay=0;];
If[aown==2&&sown==1,pay=r;];
If[aown==1&&aother==2&&sown==2,pay=pi1+c;];
If[aown==1&&aother==2&&sown==1,pay=pi1;];
If[aown==1&&aother==1&&sown==2,pay=pi2+c;];
If[aown==1&&aother==1&&sown==1,pay=pi2;];,
(* version where payoffs everything except (a1,a2) = (1,1) are normalized *)
pay=0;
If[aown==1&&aother==1,pay=theta[[s2x[{sown,sother}]]];];
If[aown==1&&sown==2,pay=pay+theta[[5]]];
];
pay
];

(* s = 2 means was inactive, a = 2 means is inactive *)
pi[a_,s_,theta_,norm_]:={
pii[a[[1]],a[[2]],s[[1]],s[[2]],theta,norm],
pii[a[[2]],a[[1]],s[[2]],s[[1]],theta,norm]
};



(* transition matrix *)
PM[sigma_]:=Map[{(1-#[[1]])*(1-#[[2]]),(1-#[[1]])*#[[2]],#[[1]]*(1-#[[2]]),#[[1]]*#[[2]]}&,sigma];
(* full strategy matrix*)
PC[sigma_]:=Module[{sigmaFull},
sigmaFull=Map[{#,1-#}&,sigma,{2}];
Table[sigmaFull[[x]][[1]][[a1]]*sigmaFull[[x]][[2]][[a2]],{x,1,4},{a1,1,2},{a2,1,2}]
];

Equilibrium solver

(* need code to solve for symmetric eq'm of game given primitives, starting value for strategies *)

(* expected error term given choice probability *)
EEps[p_]:=NIntegrate[x*PDF[NormalDistribution[0,Sqrt[var]],x],{x,InverseCDF[NormalDistribution[0,Sqrt[var]],1-p],∞}];

(* conditional value *)
Clear[EACV];
EACV[V_,theta_,norm_]:=Table[
pi[{a1,a2},x2s[x],theta,norm]+beta*V[[s2x[trans[{a1,a2},x2s[x]]]]],
{x,1,4},{a1,1,2},{a2,1,2}];

(* ex ante value function *)
Clear[EV]
EV[sigma_?NumericArrayQ,V_,theta_,norm_]:=Module[{eacv,pc},
eacv=EACV[V,theta,norm];
pc=PC[sigma];
Plus@@@(Plus@@@(pc*eacv))+Map[EEps[#]&,sigma,{2}]
];

(* symbolic value function variables *)
vSymb=Table[ToExpression["v"<>ToString[x]<>ToString[i]],{x,1,4},{i,1,2}];
vVec=Flatten[vSymb];

(* solve value function symbolically *)
VFSolve[sigma_,theta_,norm_]:=Module[{vIt,freq,bellmanEq,repl},
(*{vIt,freq}=EV[sigma,vSymb,eSim];*)
vIt=EV[sigma,vSymb,theta,norm];
bellmanEq=Flatten[vSymb==vIt];
repl=First[Solve[bellmanEq,vVec]];
N[vSymb/.repl]
];

Clear[BR];
(* best response function *)
BR[sigma_?NumericArrayQ,V_,theta_,norm_]:=Module[{sigma1,sigma2,eacv1,eacv2,eav1,eav2,deav1,deav2,br1,br2},

sigma1=Map[#[[1]]&,sigma];
sigma2=Map[#[[2]]&,sigma];

eacv1=Map[#[[1]]&,EACV[V,theta,norm],{3}];
sigma2=Map[{#,1-#}&,sigma2];
eacv1=Transpose[eacv1,{1,3,2}];
eav1=Plus@@@(sigma2*eacv1);
deav1=Map[#[[1]]-#[[2]]&,eav1];
br1=Map[CDF[NormalDistribution[0,Sqrt[var]],#]&,deav1];

eacv2=Map[#[[2]]&,EACV[V,theta,norm],{3}];
sigma1=Map[{#,1-#}&,sigma1];
eav2=Plus@@@(sigma1*eacv2);
deav2=Map[#[[1]]-#[[2]]&,eav2];
br2=Map[CDF[NormalDistribution[0,Sqrt[var]],#]&,deav2];

Table[{br1[[x]],br2[[x]]},{x,1,4}]
];
BR[sigma_?NumericArrayQ,theta_,norm_]:=BR[sigma,VFSolve[sigma,theta,norm],theta,norm];

FOC[sigma_?NumericArrayQ,theta_,norm_]:=Module[{sigma1,sigma2,sigma1full,sigma2full,V,eacv1,eacv2,eav1,eav2,dpids1,dpids2},

sigma1=Map[#[[1]]&,sigma];
sigma2=Map[#[[2]]&,sigma];
sigma1full=Map[{#,1-#}&,sigma1];
sigma2full=Map[{#,1-#}&,sigma2];

V=VFSolve[sigma,theta,norm];
Vtemp=V;

eacv1=Map[#[[1]]&,EACV[V,theta,norm],{3}];
eacv1=Transpose[eacv1,{1,3,2}];
eav1=Plus@@@(sigma2full*eacv1);
dpids1=Map[#[[1]]-#[[2]]&,eav1]-Map[InverseCDF[NormalDistribution[0,Sqrt[var]],#]&,sigma1];

eacv2=Map[#[[2]]&,EACV[V,theta,norm],{3}];
eav2=Plus@@@(sigma1full*eacv2);
dpids2=Map[#[[1]]-#[[2]]&,eav2]-Map[InverseCDF[NormalDistribution[0,Sqrt[var]],#]&,sigma2];

Table[{dpids1[[x]],dpids2[[x]]},{x,1,4}]
];


SolveSymmetric[theta_,norm_,updates_]:=Module[{iter},

sse[sigma_?NumericArrayQ]:=Module[{objval},
objval=Plus@@(Flatten[FOC[sigma,theta,norm]]^2);
If[Mod[iter,updates]==0,Print[{sigma,objval,iter}]];
objval
];

(* strategy parameters *)
sSymb={{s00,s00},{s01,s10},{s10,s01},{s11,s11}};
sVars={s00,s10,s01,s11};
sStart={{s00,0.58},{s01,.3},{s10,.84},{s11,.59}};
(* contstraint on probabilities *)
sCons=Flatten[Table[0<=sSymb[[x]][[i]]≤1,{x,1,4},{i,1,2}]];

iter=0;
Timing[
eqmB=FindMinimum[{sse[sSymb],And@@sCons},sStart,MaxIterations->25,AccuracyGoal->4,PrecisionGoal->4,EvaluationMonitor:>iter++];
];
eqmB
];

Solve for Equilibrium

(* solve for PSD's symmetric equilibrium *)
eqmBase=SolveSymmetric[basetheta,1,50]
VBase=Vtemp;

(* solve for equilibrium of counterfactual game *)
dpi={0,-.25,0,0};
cftheta=basetheta+dpi;
eqmCF=SolveSymmetric[cftheta,1,50]

eqmBase
eqmCF

(* joint CCPs in observed equilibrium *)
sigmaBase=sSymb/.eqmBase[[2]];
sigmaCF=sSymb/.eqmCF[[2]];

VBase
Vtemp

Recover primitives

(* Hotz-Miller with normal errors, varia;nce=1 *)
Clear[PhiNormal];
PhiNormal[p_?NumericQ]:=InverseCDF[NormalDistribution[0,1],p];

(* Arcidiacono-Miller with normal errors, variance=1 *)
Clear[PsiNormal]
PsiNormal[p_?NumericQ,a_]:=If[a==1,
NIntegrate[Max[eps1,-PhiNormal[p]]*PDF[NormalDistribution[0,1],eps1],{eps1,-∞,∞}],
NIntegrate[Max[PhiNormal[p]+eps1,0]*PDF[NormalDistribution[0,1],eps1],{eps1,-∞,∞}]
];

(* expected error term added to action 1*)
EepsNormal[p_]:=NIntegrate[eps1*PDF[NormalDistribution[0,1],eps1],{eps1,-PhiNormal[p],∞}]/p;

(* convert payoff vector to array indexed by own action, other action, own state, other state *)
Vec2Array[vec_]:=Partition[Partition[Partition[vec,2],2],2];

(* convert payoff array to vector *)
Array2Vec[array_]:=Flatten[array];

piSymbArray=Table[ToExpression["pi"<>ToString[a1]<>ToString[a2]<>ToString[s1]<>ToString[s2]],{a1,1,2},{a2,1,2},{s1,1,2},{s2,1,2}];
piSymbVec=Array2Vec[piSymbArray];
(* true values *)
truePi=Flatten[Table[piSymbArray[[a1]][[a2]][[s1]][[s2]]->pii[a1,a2,s1,s2,basetheta,1],{a1,1,2},{a2,1,2},{s1,1,2},{s2,1,2}]];
Vec2Array[piSymbVec]==piSymbArray

(* transition matrix, integrated over actions *)
PM[sigma_]:=Table[
{(1-sigma[[x]][[1]])*(1-sigma[[x]][[2]]),
(1-sigma[[x]][[1]])*sigma[[x]][[2]],
sigma[[x]][[1]]*(1-sigma[[x]][[2]]),
sigma[[x]][[1]]*sigma[[x]][[2]]
},{x,1,4}];

(* transition matrix, conditional on first action *)
PM1[sigma_]:=Map[{0,0,(1-#[[2]]),#[[2]]}&,sigma];

(* transition matrix, conditional on second action *)
PM2[sigma_]:=Map[{(1-#[[2]]),#[[2]],0,0}&,sigma];

Epi[sigma_]:=Module[{piSymbArray,expectedFlowProfits,expectedFlowProfits1,expectedFlowProfits2,expectedShocks,Vsymb,cvDiff,transMat,transMat1,transMat2,stratMat,ownCCP,ccpOther},

(* transition matrices *)
transMat=PM[sigma];
transMat1=PM1[sigma];
transMat2=PM2[sigma];

(* probability of pairs of actions *)
stratMat=Table[
If[a1==1,sigma[[s2x[{s1,s2}]]][[1]],1-sigma[[s2x[{s1,s2}]]][[1]]]*
If[a2==1,sigma[[s2x[{s1,s2}]]][[2]],1-sigma[[s2x[{s1,s2}]]][[2]]],
{a1,1,2},{a2,1,2},{s1,1,2},{s2,1,2}];

(* other player's strategy *)
ccpOther=
Table[
If[a2==1,sigma[[s2x[{s1,s2}]]][[2]],1-sigma[[s2x[{s1,s2}]]][[2]]],
{a2,1,2},{s1,1,2},{s2,1,2}];

piSymbArray=Table[ToExpression["pi"<>ToString[a1]<>ToString[a2]<>ToString[s1]<>ToString[s2]],{a1,1,2},{a2,1,2},{s1,1,2},{s2,1,2}];

(* expected flow profits for each state *)
expectedFlowProfits=Plus@@(Plus@@(piSymbArray*stratMat));
(* expected flow profits, with a=1 *)
expectedFlowProfits1=Plus@@(piSymbArray[[1]]*ccpOther);
(* expected flow profits, with a=2 *)
expectedFlowProfits2=Plus@@(piSymbArray[[2]]*ccpOther);

(* convert from (s1,s2) to x indices  *)
expectedFlowProfits=Flatten[Table[Extract[expectedFlowProfits,{x2s[x]}],{x,1,4}]];
expectedFlowProfits1=Flatten[Table[Extract[expectedFlowProfits1,{x2s[x]}],{x,1,4}]];
expectedFlowProfits2=Flatten[Table[Extract[expectedFlowProfits2,{x2s[x]}],{x,1,4}]];

{expectedFlowProfits,expectedFlowProfits1,expectedFlowProfits2}
];

CVDiff[sigma_]:=Module[{piSymbArray,expectedFlowProfits,expectedFlowProfits1,expectedFlowProfits2,expectedShocks,Vsymb,cvDiff,transMat,transMat1,transMat2,stratMat,ownCCP,ccpOther},

(* transition matrices *)
transMat=PM[sigma];
transMat1=PM1[sigma];
transMat2=PM2[sigma];

(* individual player's CCP's in observed equilibrium *)
ownCCP=Map[{#[[1]],1-#[[1]]}&,sigma];

(* expected error terms for each state *)
(*expectedLogit=Plus@@@(ownCCP*(EulerGamma-Log[sigmaBase]));*)
(* note: shock only to first action *)
expectedShocks=Plus@@@(ownCCP*Transpose[{Map[EepsNormal[#]&,Transpose[sigma][[1]]],Table[0,{4}]}]);

(* flow profits *)
{expectedFlowProfits,expectedFlowProfits1,expectedFlowProfits2}=Epi[sigma];

(* ex-ante value function in terms of utility function *)
Vsymb=Inverse[IdentityMatrix[4]-beta*transMat].(expectedFlowProfits+expectedShocks);

(* conditional values difference *)
expectedFlowProfits1-expectedFlowProfits2+beta*(transMat1-transMat2).Vsymb
];

(* HM inversion *) 
ownCCPBase=Map[{#[[1]],1-#[[1]]}&,sigmaBase];
(*eqConditions=Map[Log[#[[1]]/#[[2]]]&,ownCCP]-cvDiff;*)
Map[PhiNormal[#[[1]]]&,ownCCPBase]-(CVDiff[sigmaBase]/.truePi)

ObjFn[theta_,sigma_,norm_]:=Module[{piRep,cvDiff},
cvDiff=CVDiff[sigma];
piRep=Flatten[Table[piSymbArray[[a1]][[a2]][[s1]][[s2]]->pii[a1,a2,s1,s2,theta,norm],{a1,1,2},{a2,1,2},{s1,1,2},{s2,1,2}]];
errs=Map[PhiNormal[#[[1]]]&,ownCCPBase]-(cvDiff/.piRep);
Plus@@Map[#^2&,errs]
];

thetaSymb={csymb,pi1symb,pi2symb};
startingVals={-.15,1.1,-1.1};
cons={pi1symb≤100,pi1symb≥-100,pi2symb≤100,pi2symb≥-100};
starter=Table[{thetaSymb[[i]],startingVals[[i]]},{i,1,Length[thetaSymb]}];
FindMinimum[{ObjFn[{0,csymb,pi1symb,pi2symb},sigmaBase,1],cons},starter]

thetaSymb={arg1,arg2,arg3,arg4};
startingVals={-.15,-.5,2,1};
starter=Table[{thetaSymb[[i]],startingVals[[i]]},{i,1,Length[thetaSymb]}];
wrongrep=FindMinimum[ObjFn[{arg1,arg2,arg3,arg4,0},sigmaBase,2],starter]

alttheta={arg1,arg2,arg3,arg4,0}/.wrongrep[[2]];
sigmaWrongBase=SolveSymmetric[alttheta,2,50];
sigmaWrongCF=SolveSymmetric[alttheta+{0,0,0,0,dpi[[2]]},2,50];

eqmBase
eqmCF
sigmaWrongBase
sigmaWrongCF

(* have a look at payoffs for alternative normalization *)
Table[pi[{a1,a2},{s1,s2},alttheta,2],{a1,1,2},{a2,1,2},{s1,1,2},{s2,1,2}]
Table[pi[{a1,a2},{s1,s2},basetheta,1],{a1,1,2},{a2,1,2},{s1,1,2},{s2,1,2}]

(* make sure cf makes the same difference in payoffs for both normalizations *)
Table[pi[{a1,a2},{s1,s2},basetheta,1]-pi[{a1,a2},{s1,s2},basetheta+dpi,1],{a1,1,2},{a2,1,2},{s1,1,2},{s2,1,2}]-Table[pi[{a1,a2},{s1,s2},alttheta,2]-pi[{a1,a2},{s1,s2},alttheta+{0,0,0,0,dpi[[2]]},2],{a1,1,2},{a2,1,2},{s1,1,2},{s2,1,2}]

(* are the counterfactual equilibria the same? *)
(sSymb/.eqmCF[[2]])-(sSymb/.sigmaWrongCF[[2]])

Table formatting functions

NormalNumber[x_]:=ToString[x,InputForm];
tableFormat[x_,dig_]:=Module[{str},
Switch[x,
_?NumberQ,
If[
Round[x]==x,
NormalNumber[Round@x],
str=ToString[NumberForm[x,dig]];
If[Last[Characters[str]]==".",
StringJoin@Drop[Characters[str],-1],
str
]
],
_,ToString[x]
]
];

FormatNumbers[inMat_,format_]:=Map[
Switch[#,
_?NumberQ,Integerise[NumberForm[#,format]],
Null,"",
_,ToString[#]
]&,
inMat,{2}];

SmallRow[row_]:=Map["\\small{"<>#<>"}"&,row];
ItalicRow[row_]:=Map[If[#!="","\\emph{"<>#<>"}",""]&,row];
AddParens[row_?VectorQ]:=Map[If[#!="","("<>ToString[#]<>")",""]&,row];

AddParens[num_?NumberQ]:="("<>ToString[num]<>")";
AddParens[num_?StringQ]:="("<>num<>")";

AddParens[inMat_,rows_]:=Module[{i,row,outMat},
outMat=inMat;
For[i=1,i<=Length[rows],i++,
row=rows[[i]];
outMat=ReplacePart[outMat,row->
Map[If[#!="","("<>ToString[#]<>")",""]&,outMat[[row]]]
];
];
outMat
];

AddCol[inMat_,posList_?ArrayQ]:=Transpose[Insert[Transpose[inMat],Table[" ",{Length[inMat]}],posList]];
AddCol[inMat_,col_,posList_?ArrayQ]:=Transpose[Insert[Transpose[inMat],col,posList]];
AddCol[inMat_,pos_]:=Transpose[Insert[Transpose[inMat],Table[" ",{Length[inMat]}],pos]];
AddCol[inMat_,col_,pos_]:=Transpose[Insert[Transpose[inMat],col,pos]];

AddRow[inMat_,posList_?ArrayQ]:=Insert[inMat,Table[" ",{Length[inMat[[1]]]}],posList];
AddRow[inMat_,row_,posList_?ArrayQ]:=Insert[inMat,row,posList];
AddRow[inMat_,pos_]:=Insert[inMat,Table[" ",{Length[inMat[[1]]]}],pos];
AddRow[inMat_,row_,pos_]:=Insert[inMat,row,pos];

ToTex[xx_ ^yy_]:="$"<>ToTex[xx]<>"^{"<>ToTex[yy]<>"}$";
ToTex[Subscript[xx_, yy_] ]:="$"<>ToTex[xx]<>"_{"<>ToTex[yy]<>"}$";
ToTex[xx_?(#==α||#=="α"&)]:="$\\alpha$";
ToTex[xx_?(#==β||#=="β"&)]:="$\\beta$";
ToTex[xx_?(#==γ||#=="γ"&)]:="$\\gamma$";
ToTex[xx_?(#==δ||#=="δ"&)]:="$\\delta$";
ToTex[xx_?(#==ϵ||#=="ϵ"&)]:="$\\varepsilon$";
ToTex[xx_?(#==ζ||#=="ζ"&)]:="$\\zeta$";
ToTex[xx_?(#==η||#=="η"&)]:="$\\eta$";
ToTex[xx_?(#==θ||#=="θ"&)]:="$\\theta$";
ToTex[xx_?(#==κ||#=="κ"&)]:="$\\kappa$";
ToTex[xx_?(#==λ||#=="λ"&)]:="$\\lambda$";
ToTex[xx_?(#==μ||#=="μ"&)]:="$\\mu$";
ToTex[xx_?(#==ν||#=="ν"&)]:="$\\nu$";
ToTex[xx_?(#==ξ||#=="ξ"&)]:="$\\xi$";
ToTex[xx_?(#==π||#=="π"&)]:="$\\pi$";
ToTex[xx_?(#==ρ||#=="ρ"&)]:="$\\rho$";
ToTex[xx_?(#==σ||#=="σ"&)]:="$\\sigma$";
ToTex[xx_?(#==τ||#=="τ"&)]:="$\\tau$";
ToTex[xx_?(#==ϕ||#=="ϕ"&)]:="$\\phi$";
ToTex[xx_?(#==φ||#=="φ"&)]:="$\\varphi$";
ToTex[xx_?(#==χ||#=="χ"&)]:="$\\xi$";
ToTex[xx_?(#==ψ||#=="ψ"&)]:="$\\psi$";
ToTex[xx_?(#==ω||#=="ω"&)]:="$\\omega$";
ToTex[Null]:="";
ToTex[xx_?NumberQ]:=ToString[xx];
ToTex[xx_?AtomQ]:=ToString[xx];
ToTex[xx_]:=ToString[xx];

TexRow[list_]:=StringJoin[Table[ToTex[list[[i]]]<>" & ",
{i,1,Length[list]-1}]<>ToTex[list[[Length[list]]]]<>" \\\\ \n"];

Table2Tex[mat_]:=StringJoin[TexRow/@mat]

Tables

(* table of CCPs *)
pTableData=Transpose[{
Transpose[sSymb][[1]]/.eqmBase[[2]],
Transpose[BR[sigmaBase,basetheta+dpi,1]][[1]],
Transpose[BR[sigmaBase,alttheta+{0,0,0,0,dpi[[2]]},2]][[1]],
Transpose[sSymb][[1]]/.eqmCF[[2]],
Transpose[sSymb][[1]]/.sigmaWrongCF[[2]]
}];
pTableData=Map[tableFormat[#,3]&,pTableData,{2}];
pTableData//MatrixForm
collabels1={"\\multicolumn{2}{c}{$P \\left( \\mbox{active} | s_{i},s_{i}\\right)$}","Baseline","\\multicolumn{2}{c}{CF -- fixed opponent}"," " ,"\\multicolumn{2}{c}{CF -- equilibrium}" };
collabels2={" "," "," ","Model 1","Model 2"," " ,"Model 1","Model 2" };
collabels3={"$s_{i}$","$s_{-i}$"," "," "," "," "," " ," " };
rowlabels1=ToString/@{0,0,1,1};
rowlabels2=ToString/@{0,1,0,1};
pTable=Transpose[Join[{rowlabels1,rowlabels2},Transpose[pTableData]]];
pTable=Map[Insert[#," ",6]&,pTable];
pTable=
"\\begin{tabular}{cccccccc}\n"<>
"\\hline \\hline \n"<>
TexRow[collabels1]<>"\\cline{4-5} \\cline{7-8} "<>
TexRow[collabels2]<>"\\hline "<>
TexRow[collabels3]<>"\\cline{1-2} "<>
Table2Tex[pTable]<>" \\hline \\\\ \n"<> 
" \\end{tabular} "
Export["F:/Dropbox/land_values/IJIO/pTable.tex",pTable,"Text"]

(* table of payoffs *)
uTableData=Table[
Flatten[
Join[Table[pii[a1,a2,s1,s2,basetheta,1],{a1,2,1,-1},{a2,2,1,-1}],
Table[pii[a1,a2,s1,s2,alttheta,2],{a1,2,1,-1},{a2,2,1,-1}]
]
],
{s1,2,1,-1},{s2,2,1,-1}];
uTableData=Flatten[uTableData,1];
uTableData=Map[tableFormat[#,3]&,uTableData,{2}];
uTableData//MatrixForm
collabels1={"\\multicolumn{2}{c}{$\\pi \\left(a_{i},a_{i},s_{i},s_{-i}\\right)$}","\\multicolumn{4}{c}{Model 1}"," ","\\multicolumn{4}{c}{Model 2}"};
collabels2={" "," ","\\multicolumn{4}{c}{$\\left(a_{i},a_{-i}\\right)=$}"," ","\\multicolumn{4}{c}{$\\left(a_{i},a_{-i}\\right)=$}"};
collabels3={"$s_{i}$","$s_{-i}$","$\\left(0,0\\right)$","$\\left(0,1\\right)$","$\\left(1,0\\right)$","$\\left(1,1\\right)$"," ","$\\left(0,0\\right)$","$\\left(0,1\\right)$","$\\left(1,0\\right)$","$\\left(1,1\\right)$"};
rowlabels1=ToString/@{0,0,1,1};
rowlabels2=ToString/@{0,1,0,1};
uTable=Transpose[Join[{rowlabels1,rowlabels2},Transpose[uTableData]]];
uTable=Map[Insert[#," ",7]&,uTable];
uTable=
"\\begin{tabular}{ccccccccccc}\n"<>
"\\hline \\hline \n"<>
TexRow[collabels1]<>"\\cline{3-6} \\cline{8-11} "<>
TexRow[collabels2]<>
TexRow[collabels3]<>"\\hline "<>
Table2Tex[uTable]<>" \\hline \\\\ \n"<> 
" \\end{tabular} "
Export["F:/Dropbox/land_values/IJIO/uTable.tex",uTable,"Text"]