#====================================================================#
# This is the main R code for the simulation part of the paper:
# "A Distributional Synthetic Control Method for Policy Evaluation"
# This version: January 9, 2020
#====================================================================#
rm(list=ls())
library(Matrix)
library(quadprog)

#==========================#
# Part I. Define functions # 
#==========================#

########################################################
# Function: Synthetic Control by Quadratic Programming #
########################################################
# Yseq: full sequence of treatment unit
# Xseq: full sequences of potential control units
SynthCtl=function(Yin,Xin,Yz,Xz){  
# Quadratic programming for minimizing ||Y-X*w||2
Yin=as.matrix(Yin)
Xin=as.matrix(Xin)
Yz=as.matrix(Yz)
Xz=as.matrix(Xz)
Y=c(Yz,Yin)
X=rbind(Xz,Xin)
S=length(Y)             # length of data
n=ncol(X)               # number of potential control units
wMtx=matrix(0,n,1)      # synthetic weights
Synth=matrix(0,S,1)     # Synthetic-control sequence
Dmat=t(X)%*%X/S
dvec=t(X)%*%Y/S
Amat=rbind(matrix(1,1,n),diag(n)) 
Amat=t(Amat)
bvec0=c(1,matrix(0,n,1))
pd_Dmat=nearPD(Dmat)    # Modification: near positive definite matrix
qp=solve.QP(Dmat=as.matrix(pd_Dmat$mat),dvec,Amat,bvec=bvec0,meq=1) 
#qp=solve.QP(Dmat,dvec,Amat,bvec=bvec0,meq=1) 
wMtx=qp$solution           
Synth=Xin%*%wMtx
MSE=mean((Yin-Synth)^2)
list(res1=wMtx,res2=Synth,res3=MSE)}

#####################################
# Function: Data generating process #
#####################################
# Inputs:
# N: number of units (one treated unit + (N-1) control units) 
# n: number of subuits for each unit
# T: length of time series
# T0: end of the pre-intervetion period
# Alpha: coefficients of observed factors
# beta: coefficients of structural breaks
# tau: a vector probability-index of quantiles
# Outputs:
# Y: subunit-level panel data for each unit
# Q: quantile time series for each unit
# x: predictors for each unit
DGP=function(N,n,T,T0,Alpha,xi,tau){  
m=length(tau)
g=rlnorm(T)
x=rlnorm(N)
alphai0=matrix(0,T,N)  # constant
alphai1=matrix(0,T,N)  # time trend coefficient
gamma=rlnorm(N)
delta=matrix(0,T,1)
policy=matrix(0,T,1)
t=1
while(t<=T){
alphai0[t,]=Alpha[,1]	
if(t<=T0){alphai1[t,]=Alpha[,2]}
if(t>T0){alphai1[t,]=Alpha[,2]+xi[t,]; delta[t]=2*(t-T0)}
t=t+1}
Y=matrix(0,T,(N*n))
Q=matrix(0,T,(N*m))
i=1
while(i<=N){
#if(i==1){alphai=as.matrix((1+alpha)*c(1,0.5))}	
#if(i>1){alphai=as.matrix(c(1,0.5))}	
if(i==1){eff=delta}	
if(i>1){eff=matrix(0,T,1)}	
err=rnorm(T)
t=1
while(t<=T){
tauY=runif(n)
zit=rlnorm(1)
zit=1
Y[t,((i-1)*n+1):(i*n)]=(g[t]+alphai0[t,i]+alphai1[t,i]*t+x[i]+gamma[i]*zit+eff[t])*exp(tauY)+err[t]
Q[t,((i-1)*m+1):(i*m)]=quantile(Y[t,((i-1)*n+1):(i*n)],tau)
t=t+1}
i=i+1}
list(res1=Y,res2=Q,res=x)}

#=================================#
# Part II. Monte Carlo simulation # 
#=================================#
dirpath="E:/Research/SyntheticControl/Revision2/Simulation/"

# Probability-index vector of quantiles
tau=c(0.01,0.05,0.1,0.3,0.5,0.7,0.9,0.95,0.99)
m=length(tau)


# number of replications
R=1000

T_type=1
while(T_type<=2){
N_type=1
while(N_type<=2){
n_type=2
while(n_type<=2){
M_type=1
while(M_type<=1){
S_type=2
while(S_type<=2){

# Sample sizes
if(T_type==1){T0=50}  # shorter pre-intervention duration
if(T_type==2){T0=100} # longer pre-intervention duration
if(N_type==1){N=50}   # smaller number of units
if(N_type==2){N=100}  # larger number of units
if(n_type==1){n=50}   # smaller number of subunits
if(n_type==2){n=100}  # larger number of subunits
T=T0+10               # length of time series


# Poor matching 
alphaC=seq(1/N,(1-1/N),1/N)
if(M_type==1){alpha1=0.5} # no poor matching (the treated unit is in the convex hull)
if(M_type==2){alpha1=1.5} # poor matching (the treated unit is out of the convex hull)
alpha=c(alpha1,alphaC)
Alpha=cbind(alpha,alpha/2)

# Structural break
if(S_type==1){      # no structural break
xi=matrix(0,T,N)} 
if(S_type==2){      # structural break
xi=matrix(0,T,N)
t=1
while(t<=T){
if(t>T0){xi[t,]=0.1*Alpha[,2]*(t-T0)}
t=t+1}} 


# Matrices to be saved
Mtrx_Trt=matrix(0,T,(R*m))   # Simulated acutal quantial time series
Mtrx_Syn0=matrix(0,T,(R*m))  # Simulated counterfactutal quantial time series (traditional)
Mtrx_Syn=matrix(0,T,(R*m))   # Simulated counterfactutal quantial time series (proposed)

# random seed
set.seed(123)  
r=1
while(r<=R){

#-----------------------------------------------------------------------------------	
# (1) Data generating process of the subunit-level panels
#-----------------------------------------------------------------------------------	
res=DGP(N,n,T,T0,Alpha,xi,tau)
Y=res[[1]]
Q=res[[2]]
XData=res[[3]]

## The subunit time series and the quantile time series of the treated unit
#i=1
#par(mfrow=c(1,2),mar=c(2,2,2,2))
#matplot(Y[,(n*(i-1)+1):(n*i)],type="l")
#abline(v=T0)
#matplot(Q[,(m*(i-1)+1):(m*i)],type="l")
#abline(v=T0)
#par(mfrow=c(1,1))
## The subunit time series and the quantile time series of the potential control units
#i=2
#while(i<=N){
#par(mfrow=c(1,2),mar=c(2,2,2,2))
#matplot(Y[,(n*(i-1)+1):(n*i)],type="l")
#abline(v=T0)
#matplot(Q[,(m*(i-1)+1):(m*i)],type="l")
#abline(v=T0)
#par(mfrow=c(1,1))
#i=i+1}

#-----------------------------------------------------------------------------------	
# (2) Synthetic control analysis
#-----------------------------------------------------------------------------------	
# Define time indices
a=1      # the start of the pre-intervention period  
b=T0     # the end of the pre-intervention period
c=T0+1   # the start of the post-intervention period (policy-change point) 
d=T      # the end of the post-intervention period    

q=1
while(q<=(length(tau))){

Ymtrx=matrix(0,T,N)      # the tau-quantile panel
i=1
while(i<=N){
Ymtrx[,i]=Q[,(m*(i-1)+q)]
i=i+1}

# the treated (1st) unit
trtPre=as.matrix(Ymtrx[a:b,1])   # pre-intervention time series
trtPost=as.matrix(Ymtrx[c:d,1])  # post-intervention time series
trt_x=as.matrix(XData[1])        # predictor

# the (N-1) potential control units
pcuPre=Ymtrx[a:b,-1]             # pre-intervention time series
pcuPost=Ymtrx[c:d,-1]            # post-intervention time series
pcu_x=t(XData[-1])               # predictor

# h=(1,t) for controlling the poor-matching issue
k=1
F=matrix(0,(b-a+1),(k+1))
t=1
while(t<=(b-a+1)){
j=0
while(j<=k){
F[t,(j+1)]=t^j
j=j+1}
t=t+1}
X=F

# OLS estimation
trt_alpha=solve(t(X)%*%X)%*%(t(X)%*%trtPre)
pcu_alpha=solve(t(X)%*%X)%*%(t(X)%*%pcuPre)
trtPre_fit=X%*%trt_alpha      # OLS fitted values
pcuPre_fit=X%*%pcu_alpha
trtPre_res=trtPre-trtPre_fit  # OLS residuals
pcuPre_res=pcuPre-pcuPre_fit

# Optimal weights of the traditional synthetic control (ignores the poor matching)
RES0=SynthCtl(trtPre,pcuPre,trt_x,pcu_x)
optw0=RES0[[1]]  # optimal weights obtained from the quadratic programming 
syth0=RES0[[2]]  # the pre-intervention synthetic-control data

# Optimal weights of the proposed synthetic control (modifies the poor matching)
RES=SynthCtl(trtPre_res,pcuPre_res,trt_x,pcu_x)
optw=RES[[1]]    # optimal weights obtained from the quadratic programming 
syth=RES[[2]]    # the pre-intervention synthetic-control data

# Post-intervention synthetic controls
synPost0=matrix(0,length(trtPost),1)               
synPost=matrix(0,length(trtPost),1)                
FPost=matrix(0,(d-c+1),(k+1))
s=1
while(s<=(d-c+1)){
j=0
while(j<=k){
FPost[s,(j+1)]=(c-a+s)^j
j=j+1}
s=s+1}
XPost=FPost
trtPost_fit=XPost%*%trt_alpha    # OLS fitted values
pcuPost_fit=XPost%*%pcu_alpha
i=1
while(i<=(N-1)){
synPost0=synPost0+pcuPost[,i]*optw0[i]
synPost=synPost+(pcuPost[,i]-pcuPost_fit[,i])*optw[i]
i=i+1}
synPre=syth+trtPre_fit
synPost=synPost+trtPost_fit
trt_fit=c(trtPre_fit,trtPost_fit)

#-----------------------------------------------------------------------------------	
# (3) Acutal and counterfactual quantile time series of the treated unit
#-----------------------------------------------------------------------------------	
Trt=c(trtPre,trtPost)        # Acutal quantial time series
MySyn0=c(syth0,synPost0)     # Counterfactual quantial time series (traditional)
MySyn=c(synPre,synPost)      # Counterfactual quantial time series (proposed)

mt=paste0("S=",S_type,"M=",M_type,"n=",n_type,"N=",N_type,"T=",T_type,"r=",r,"q=",q)                                          
plot(Trt,main=mt,type="l",lty=1,lwd=4,col=2,cex=1,xlab=" ",ylab=" ")
lines(MySyn0,type="l",lwd=2,lty=2,col=4)
lines(MySyn,type="l",lwd=2,lty=1,col=1)
abline(v=T0)
legend("topleft",legend=c("Treated","Traditional","Proposed"),lty=c(1,2,1),lwd=c(4,2,2),col=c(2,4,1),bty="n",cex=1,y.intersp=1)

Mtrx_Trt[,((r-1)*m+q)]=Trt
Mtrx_Syn0[,((r-1)*m+q)]=MySyn0
Mtrx_Syn[,((r-1)*m+q)]=MySyn

q=q+1}
r=r+1}

if(S_type==1 & M_type==1 & T_type==1 & N_type==1 & n_type==1){write.csv(Mtrx_Trt, file=paste0(dirpath,"Trt_S1M1T1N1n1.csv")); write.csv(Mtrx_Syn0, file=paste0(dirpath,"Syn0_S1M1T1N1n1.csv")); write.csv(Mtrx_Syn, file=paste0(dirpath,"Syn_S1M1T1N1n1.csv")); } 
if(S_type==1 & M_type==1 & T_type==1 & N_type==1 & n_type==2){write.csv(Mtrx_Trt, file=paste0(dirpath,"Trt_S1M1T1N1n2.csv")); write.csv(Mtrx_Syn0, file=paste0(dirpath,"Syn0_S1M1T1N1n2.csv")); write.csv(Mtrx_Syn, file=paste0(dirpath,"Syn_S1M1T1N1n2.csv")); } 
if(S_type==1 & M_type==1 & T_type==1 & N_type==2 & n_type==1){write.csv(Mtrx_Trt, file=paste0(dirpath,"Trt_S1M1T1N2n1.csv")); write.csv(Mtrx_Syn0, file=paste0(dirpath,"Syn0_S1M1T1N2n1.csv")); write.csv(Mtrx_Syn, file=paste0(dirpath,"Syn_S1M1T1N2n1.csv")); } 
if(S_type==1 & M_type==1 & T_type==1 & N_type==2 & n_type==2){write.csv(Mtrx_Trt, file=paste0(dirpath,"Trt_S1M1T1N2n2.csv")); write.csv(Mtrx_Syn0, file=paste0(dirpath,"Syn0_S1M1T1N2n2.csv")); write.csv(Mtrx_Syn, file=paste0(dirpath,"Syn_S1M1T1N2n2.csv")); } 
if(S_type==1 & M_type==1 & T_type==2 & N_type==1 & n_type==1){write.csv(Mtrx_Trt, file=paste0(dirpath,"Trt_S1M1T2N1n1.csv")); write.csv(Mtrx_Syn0, file=paste0(dirpath,"Syn0_S1M1T2N1n1.csv")); write.csv(Mtrx_Syn, file=paste0(dirpath,"Syn_S1M1T2N1n1.csv")); } 
if(S_type==1 & M_type==1 & T_type==2 & N_type==1 & n_type==2){write.csv(Mtrx_Trt, file=paste0(dirpath,"Trt_S1M1T2N1n2.csv")); write.csv(Mtrx_Syn0, file=paste0(dirpath,"Syn0_S1M1T2N1n2.csv")); write.csv(Mtrx_Syn, file=paste0(dirpath,"Syn_S1M1T2N1n2.csv")); } 
if(S_type==1 & M_type==1 & T_type==2 & N_type==2 & n_type==1){write.csv(Mtrx_Trt, file=paste0(dirpath,"Trt_S1M1T2N2n1.csv")); write.csv(Mtrx_Syn0, file=paste0(dirpath,"Syn0_S1M1T2N2n1.csv")); write.csv(Mtrx_Syn, file=paste0(dirpath,"Syn_S1M1T2N2n1.csv")); } 
if(S_type==1 & M_type==1 & T_type==2 & N_type==2 & n_type==2){write.csv(Mtrx_Trt, file=paste0(dirpath,"Trt_S1M1T2N2n2.csv")); write.csv(Mtrx_Syn0, file=paste0(dirpath,"Syn0_S1M1T2N2n2.csv")); write.csv(Mtrx_Syn, file=paste0(dirpath,"Syn_S1M1T2N2n2.csv")); } 
if(S_type==1 & M_type==2 & T_type==1 & N_type==1 & n_type==1){write.csv(Mtrx_Trt, file=paste0(dirpath,"Trt_S1M2T1N1n1.csv")); write.csv(Mtrx_Syn0, file=paste0(dirpath,"Syn0_S1M2T1N1n1.csv")); write.csv(Mtrx_Syn, file=paste0(dirpath,"Syn_S1M2T1N1n1.csv")); } 
if(S_type==1 & M_type==2 & T_type==1 & N_type==1 & n_type==2){write.csv(Mtrx_Trt, file=paste0(dirpath,"Trt_S1M2T1N1n2.csv")); write.csv(Mtrx_Syn0, file=paste0(dirpath,"Syn0_S1M2T1N1n2.csv")); write.csv(Mtrx_Syn, file=paste0(dirpath,"Syn_S1M2T1N1n2.csv")); } 
if(S_type==1 & M_type==2 & T_type==1 & N_type==2 & n_type==1){write.csv(Mtrx_Trt, file=paste0(dirpath,"Trt_S1M2T1N2n1.csv")); write.csv(Mtrx_Syn0, file=paste0(dirpath,"Syn0_S1M2T1N2n1.csv")); write.csv(Mtrx_Syn, file=paste0(dirpath,"Syn_S1M2T1N2n1.csv")); } 
if(S_type==1 & M_type==2 & T_type==1 & N_type==2 & n_type==2){write.csv(Mtrx_Trt, file=paste0(dirpath,"Trt_S1M2T1N2n2.csv")); write.csv(Mtrx_Syn0, file=paste0(dirpath,"Syn0_S1M2T1N2n2.csv")); write.csv(Mtrx_Syn, file=paste0(dirpath,"Syn_S1M2T1N2n2.csv")); } 
if(S_type==1 & M_type==2 & T_type==2 & N_type==1 & n_type==1){write.csv(Mtrx_Trt, file=paste0(dirpath,"Trt_S1M2T2N1n1.csv")); write.csv(Mtrx_Syn0, file=paste0(dirpath,"Syn0_S1M2T2N1n1.csv")); write.csv(Mtrx_Syn, file=paste0(dirpath,"Syn_S1M2T2N1n1.csv")); } 
if(S_type==1 & M_type==2 & T_type==2 & N_type==1 & n_type==2){write.csv(Mtrx_Trt, file=paste0(dirpath,"Trt_S1M2T2N1n2.csv")); write.csv(Mtrx_Syn0, file=paste0(dirpath,"Syn0_S1M2T2N1n2.csv")); write.csv(Mtrx_Syn, file=paste0(dirpath,"Syn_S1M2T2N1n2.csv")); } 
if(S_type==1 & M_type==2 & T_type==2 & N_type==2 & n_type==1){write.csv(Mtrx_Trt, file=paste0(dirpath,"Trt_S1M2T2N2n1.csv")); write.csv(Mtrx_Syn0, file=paste0(dirpath,"Syn0_S1M2T2N2n1.csv")); write.csv(Mtrx_Syn, file=paste0(dirpath,"Syn_S1M2T2N2n1.csv")); } 
if(S_type==1 & M_type==2 & T_type==2 & N_type==2 & n_type==2){write.csv(Mtrx_Trt, file=paste0(dirpath,"Trt_S1M2T2N2n2.csv")); write.csv(Mtrx_Syn0, file=paste0(dirpath,"Syn0_S1M2T2N2n2.csv")); write.csv(Mtrx_Syn, file=paste0(dirpath,"Syn_S1M2T2N2n2.csv")); } 

if(S_type==2 & M_type==1 & T_type==1 & N_type==1 & n_type==1){write.csv(Mtrx_Trt, file=paste0(dirpath,"Trt_S2M1T1N1n1.csv")); write.csv(Mtrx_Syn0, file=paste0(dirpath,"Syn0_S2M1T1N1n1.csv")); write.csv(Mtrx_Syn, file=paste0(dirpath,"Syn_S2M1T1N1n1.csv")); } 
if(S_type==2 & M_type==1 & T_type==1 & N_type==1 & n_type==2){write.csv(Mtrx_Trt, file=paste0(dirpath,"Trt_S2M1T1N1n2.csv")); write.csv(Mtrx_Syn0, file=paste0(dirpath,"Syn0_S2M1T1N1n2.csv")); write.csv(Mtrx_Syn, file=paste0(dirpath,"Syn_S2M1T1N1n2.csv")); } 
if(S_type==2 & M_type==1 & T_type==1 & N_type==2 & n_type==1){write.csv(Mtrx_Trt, file=paste0(dirpath,"Trt_S2M1T1N2n1.csv")); write.csv(Mtrx_Syn0, file=paste0(dirpath,"Syn0_S2M1T1N2n1.csv")); write.csv(Mtrx_Syn, file=paste0(dirpath,"Syn_S2M1T1N2n1.csv")); } 
if(S_type==2 & M_type==1 & T_type==1 & N_type==2 & n_type==2){write.csv(Mtrx_Trt, file=paste0(dirpath,"Trt_S2M1T1N2n2.csv")); write.csv(Mtrx_Syn0, file=paste0(dirpath,"Syn0_S2M1T1N2n2.csv")); write.csv(Mtrx_Syn, file=paste0(dirpath,"Syn_S2M1T1N2n2.csv")); } 
if(S_type==2 & M_type==1 & T_type==2 & N_type==1 & n_type==1){write.csv(Mtrx_Trt, file=paste0(dirpath,"Trt_S2M1T2N1n1.csv")); write.csv(Mtrx_Syn0, file=paste0(dirpath,"Syn0_S2M1T2N1n1.csv")); write.csv(Mtrx_Syn, file=paste0(dirpath,"Syn_S2M1T2N1n1.csv")); } 
if(S_type==2 & M_type==1 & T_type==2 & N_type==1 & n_type==2){write.csv(Mtrx_Trt, file=paste0(dirpath,"Trt_S2M1T2N1n2.csv")); write.csv(Mtrx_Syn0, file=paste0(dirpath,"Syn0_S2M1T2N1n2.csv")); write.csv(Mtrx_Syn, file=paste0(dirpath,"Syn_S2M1T2N1n2.csv")); } 
if(S_type==2 & M_type==1 & T_type==2 & N_type==2 & n_type==1){write.csv(Mtrx_Trt, file=paste0(dirpath,"Trt_S2M1T2N2n1.csv")); write.csv(Mtrx_Syn0, file=paste0(dirpath,"Syn0_S2M1T2N2n1.csv")); write.csv(Mtrx_Syn, file=paste0(dirpath,"Syn_S2M1T2N2n1.csv")); } 
if(S_type==2 & M_type==1 & T_type==2 & N_type==2 & n_type==2){write.csv(Mtrx_Trt, file=paste0(dirpath,"Trt_S2M1T2N2n2.csv")); write.csv(Mtrx_Syn0, file=paste0(dirpath,"Syn0_S2M1T2N2n2.csv")); write.csv(Mtrx_Syn, file=paste0(dirpath,"Syn_S2M1T2N2n2.csv")); } 
if(S_type==2 & M_type==2 & T_type==1 & N_type==1 & n_type==1){write.csv(Mtrx_Trt, file=paste0(dirpath,"Trt_S2M2T1N1n1.csv")); write.csv(Mtrx_Syn0, file=paste0(dirpath,"Syn0_S2M2T1N1n1.csv")); write.csv(Mtrx_Syn, file=paste0(dirpath,"Syn_S2M2T1N1n1.csv")); } 
if(S_type==2 & M_type==2 & T_type==1 & N_type==1 & n_type==2){write.csv(Mtrx_Trt, file=paste0(dirpath,"Trt_S2M2T1N1n2.csv")); write.csv(Mtrx_Syn0, file=paste0(dirpath,"Syn0_S2M2T1N1n2.csv")); write.csv(Mtrx_Syn, file=paste0(dirpath,"Syn_S2M2T1N1n2.csv")); } 
if(S_type==2 & M_type==2 & T_type==1 & N_type==2 & n_type==1){write.csv(Mtrx_Trt, file=paste0(dirpath,"Trt_S2M2T1N2n1.csv")); write.csv(Mtrx_Syn0, file=paste0(dirpath,"Syn0_S2M2T1N2n1.csv")); write.csv(Mtrx_Syn, file=paste0(dirpath,"Syn_S2M2T1N2n1.csv")); } 
if(S_type==2 & M_type==2 & T_type==1 & N_type==2 & n_type==2){write.csv(Mtrx_Trt, file=paste0(dirpath,"Trt_S2M2T1N2n2.csv")); write.csv(Mtrx_Syn0, file=paste0(dirpath,"Syn0_S2M2T1N2n2.csv")); write.csv(Mtrx_Syn, file=paste0(dirpath,"Syn_S2M2T1N2n2.csv")); } 
if(S_type==2 & M_type==2 & T_type==2 & N_type==1 & n_type==1){write.csv(Mtrx_Trt, file=paste0(dirpath,"Trt_S2M2T2N1n1.csv")); write.csv(Mtrx_Syn0, file=paste0(dirpath,"Syn0_S2M2T2N1n1.csv")); write.csv(Mtrx_Syn, file=paste0(dirpath,"Syn_S2M2T2N1n1.csv")); } 
if(S_type==2 & M_type==2 & T_type==2 & N_type==1 & n_type==2){write.csv(Mtrx_Trt, file=paste0(dirpath,"Trt_S2M2T2N1n2.csv")); write.csv(Mtrx_Syn0, file=paste0(dirpath,"Syn0_S2M2T2N1n2.csv")); write.csv(Mtrx_Syn, file=paste0(dirpath,"Syn_S2M2T2N1n2.csv")); } 
if(S_type==2 & M_type==2 & T_type==2 & N_type==2 & n_type==1){write.csv(Mtrx_Trt, file=paste0(dirpath,"Trt_S2M2T2N2n1.csv")); write.csv(Mtrx_Syn0, file=paste0(dirpath,"Syn0_S2M2T2N2n1.csv")); write.csv(Mtrx_Syn, file=paste0(dirpath,"Syn_S2M2T2N2n1.csv")); } 
if(S_type==2 & M_type==2 & T_type==2 & N_type==2 & n_type==2){write.csv(Mtrx_Trt, file=paste0(dirpath,"Trt_S2M2T2N2n2.csv")); write.csv(Mtrx_Syn0, file=paste0(dirpath,"Syn0_S2M2T2N2n2.csv")); write.csv(Mtrx_Syn, file=paste0(dirpath,"Syn_S2M2T2N2n2.csv")); } 


S_type=S_type+1}
M_type=M_type+1}
n_type=n_type+1}
N_type=N_type+1}
T_type=T_type+1}
















