function scalar csv_print_matrix (const matrix X, int N,
                                  const string fname,
                                  const strings vnames[null])
  NT = rows(X)
  k = cols(X)
  scalar T = NT / N
  scalar have_names = nelem(vnames) > 0
  string nl = sprintf("\n")
  scalar r, xrj
  
  set force_decpoint on

  outfile @fname --quiet
    printf "unit,time,"
    loop j=1..k -q
      if have_names
        printf "%s", vnames[j]
      else
        printf "v%d", j
      endif
      printf "%s", j == k ? nl : ","
    endloop
    r = 1
    loop i=1..N -q
      loop t=1..T -q
        printf "%d,%d,", i, t
        loop j=1..k -q
          xrj = X[r,j]
          if missing(xrj)
            printf "NA"
          else
            printf "%.15g", xrj
          endif
          printf "%s", j == k ? nl : ","
        endloop
        r++
      endloop
    endloop
  end outfile

  return 0
end function

function matrix pstack (int N, list L,
                        const string fname[null],
                        const strings vnames[null])

  if N < 2 || nelem(L) < 2
    funcerr "Both N and T must be at least 2"
  endif                      
                      
  if exists("vnames")
    nn = nelem(vnames)
  else
    strings vnames = array(0)
    nn = 0
  endif

  # don't omit any missing values  
  set skip_missing off
  matrix X = {L}
  set skip_missing on
  if rows(X) % N != 0
    funcerr "The number of rows must be an integer multiple of N"
  endif
  
  T = cols(X)     # number of periods
  NT = N * T      # panelized series length
  k = nelem(X)/NT # number of panelized series
  X = mshape(X', NT, k)

  if nn > 0 && nn < k
    funcerr "The vnames array is too short"
  endif

  if exists("fname")
    csv_print_matrix(X, N, fname, vnames)
  endif

  return X
end function

# quit

open byvar.csv
N = 6 # number of individuals
strings vnames = defarray("foo", "bar")
matrix X = pstack(6, dataset, "foo.csv", vnames)
open foo.csv
print -o