mata:
  real matrix t2seuc(matrix G, real scalar disttype, real scalar dropdiag, real scalar rcb) {
    // Read stata matrix into mata
//    G=st_matrix(tabmat)

    if (rows(G)!=cols(G)) {
      _error("Table isn't square")
    }
    dim = rows(G)

    if (dropdiag==1) {
      T = G :- diag(G) // only proper transitions
    } else {
      T = G
    }

    Tr = T :/ rowsum(T)
    Tc = T :/ colsum(T)

    rowtots = rowsum(T)
    coltots = colsum(T)
    rowpct = T :/ rowtots
    colpct = T :/ coltots
    // Get Pearson residuals and relative residuals
    E = (rowsum(T)*colsum(T)) :/ sum(T)
    RR = (T :- E) :/ E
    PR = (T :- E) :/ sqrt(E)


    d1 = J(dim,dim,0)
    if (rcb==1 | rcb==3) { // row or both
      for (i=1; i<=dim; i++) {
        for (j=i; j<=dim; j++) {
          if (disttype==4) { // just pij and pji
            d1[i,j] = 2 - Tr[i,j] - Tr[j,i]
          } else { // for all states
            for (k = 1; k<=dim; k++) {
              if (disttype==0) { // Percentages
                d1[i,j] = d1[i,j] + (rowpct[i,k]-rowpct[j,k])^2
              } else if (disttype==6) { // px, test a transposed version of the pct distance (poor performance) Dec 21 2023 15:05:08
                d1[i,j] = d1[i,j] + (rowpct[k,i]-rowpct[k,j])^2
              } else if (disttype==1) { // Pearson residual
                d1[i,j] = d1[i,j] + (PR[i,k]-PR[j,k])^2
              } else if (disttype==2) { // Relative residual
                d1[i,j] = d1[i,j] + (RR[i,k]-RR[j,k])^2
              } else if (disttype==3) { // chisq
                // See https://link.springer.com/referenceworkentry/10.1007/978-0-387-32833-1_53
                d1[i,j] = d1[i,j] + (rowpct[i,k]-rowpct[j,k])^2/coltots[k]
              }
            }
          }
        }
      }
    }
    d1 = d1 :+ transposeonly(d1)

    d2 = J(dim,dim,0)
    if (rcb==2 | rcb==3) { // col or both
      for (i=1; i<=dim; i++) {
        for (j=i; j<=dim; j++) {
          if (disttype==4) { // just pij and pji
            d2[i,j] = 2 - Tc[i,j] - Tc[j,i]
          } else { // for all states
            for (k = 1; k<=dim; k++) {
              if (disttype==0) { // Percentages
                d2[i,j] = d2[i,j] + (colpct[k,i]-colpct[k,j])^2
              } else if (disttype==6) { // px, test
                d2[i,j] = d2[i,j] + (colpct[i,k]-colpct[j,k])^2
              } else if (disttype==1) { // Pearson residual
                d2[i,j] = d2[i,j] + (PR[k,i]-PR[k,j])^2
              } else if (disttype==2) { // Relative residual
                d2[i,j] = d2[i,j] + (RR[k,i]-RR[k,j])^2
              } else if (disttype==3) { // chisq
                d2[i,j] = d2[i,j] + (colpct[k,i]-colpct[k,j])^2/rowtots[k]
              }
            }
          }
        }
      }
    }
    d2 = d2 :+ transposeonly(d2)

    result = sqrt(d1 :+ d2)

    if (disttype==4) {
      result = result :- diag(result) // force zeros on diagonal
    }
    return(result)
}

end

program define sdmksubs
syntax [if] [in], SUBSmat(string) [STATEvar(varname) IDvar(varname) DTYPE(string) RCB(string) RAW DROPDiag LAG(real 1) NStates(real 0)]

if ("`dropdiag'" == "") {
  local dropdiag 0
}
else {
  local dropdiag 1
}

local typeno -1
if ("`dtype'" == "pct") { // Percentages
  local typeno 0
}
if ("`dtype'" == "px") { // Percentages
  local typeno 6
}
if ("`dtype'" == "pearson") {
  local typeno 1
}
if ("`dtype'" == "relres") {
  local typeno 2
}
if ("`dtype'" == "chisq") {
  local typeno 3
}
if ("`dtype'" == "t2s") {
  local typeno 4
}
if ("`dtype'" == "flat") {
  local typeno 5
}
if ("`dtype'" == "linear") {
  local typeno 5
}

if (`typeno' == -1) {
  di in red "Type `dtype' not recognised"
  exit
  }

local rowcolboth 1 // default to row
if ("`rcb'" == "r") {
local rowcolboth 1
}
if ("`rcb'" == "c") {
local rowcolboth 2
}
if ("`rcb'" == "b") {
local rowcolboth 3
}

marksample touse

if !inlist("`dtype'", "flat", "linear") {
  if ("`idvar'"=="") {
    di in red "Option IDVAR required"
    exit
  }
  if ("`statevar'" == "") {
    di in red "State variable must be specified for type `dtype'"
    exit
  }

  tempvar rowvar seqlen
  tempname matcell

  qui by `idvar': gen `seqlen' = _N
  qui su `seqlen'
  if (`lag' > r(max) - 1) {
    di in red "Lag `lag' too large for data"
    exit
    }
  qui by `idvar': gen `rowvar'=`statevar'[_n-`lag'] if _n>1

  di "Generating substitution matrix type `dtype'"

  qui tab `rowvar' `statevar' if `touse', matcell(`matcell')

  mata: st_matrix("`subsmat'", t2seuc(st_matrix("`matcell'"), `typeno', `dropdiag', `rowcolboth'))
  if ("`dtype'" == "t2s" && `rowcolboth' != 3 ) {
    // For traditional t2s on row or col, 2-pij-pji not sqrt(2-pij-pji)
    mata: st_matrix("`subsmat'", st_matrix("`subsmat'") :^ 2)
  }
}
else {
  if ("`statevar'"=="" && `nstates'==0) {
    di in red "Either state variable or nstates() option must be specified for type `dtype'"
    exit
  }
  if (`nstates' == 0) {
    qui su `statevar'
    local dims = r(max)
  }
  else {
    local dims `nstates'
  }

  if ("`dtype'"=="flat") {
    matrix `subsmat' = J(`dims',`dims',1) - I(`dims')
  }
  else {
    matrix `subsmat' = J(`dims',`dims',0)
    forvalues x = 1/`dims' {
      forvalues y = `=`x'+1'/`dims' {
        matrix `subsmat'[`x',`y'] = abs(`x' - `y')
        matrix `subsmat'[`y',`x'] = abs(`x' - `y')
      }
    }
  }
}

if ("`raw'" == "") {
  di "Scaling to 0--1 range"
  mata: st_matrix("`subsmat'", st_matrix("`subsmat'") :/ max(st_matrix("`subsmat'")))
}

end