Riccardo "Jack" Lucchetti and Alessandro Pionati and Francesco Valentini 0.1 2026-07-28 Grouped Fixed Effects for panel data C23 pdfdoc:GFE.pdf examples chk = bal_test(y, X) if !chk funcerr "Panel is unbalanced. Not for now" endif N = $nobs/$pd # process options opts = default_opts(G) if exists(opt_in) opts = opt_in + opts endif # choose a sensible number of random starts if not given if missing(opts.randstarts) if opts.algorithm == 1 opts.randstarts = G < 8 ? 1024 : 4096 elif opts.algorithm == 2 opts.randstarts = floor(8*sqrt(G)) endif endif set stopwatch if opts.BICselect printf "%s\n", "="*64 printf " BIC model selection (Gmax = %d)\n", G printf "%s\n", "="*64 opts.BICselect = 0 verbose = opts.verbose opts.verbose = 0 mod = gfe_estimate_fixed(y, X, G, opts) bestg = G T = mod.T K = mod.nreg ssr = mod.obj s2hat = ssr / ((N - G)*T - N - K) BICmax = BICrit(ssr, N, T, G, K, s2hat) if verbose printf "%6.3g%% done\r", 100/G flush endif matrix BIC = {G, BICmax, ssr} bestBIC = BICmax loop g = (G-1) .. 2 --decr bundle b = gfe_estimate_fixed(y, X, g, opts) BICi = BICrit(b.obj, N, T, g, K, s2hat) if verbose printf "%6.3g%% done\r", 100*(G-g+1)/G flush endif BIC |= {g, BICi, b.obj} if BICi < bestBIC mod = b bestg = g bestBIC = BICi endif endloop if verbose printf " 100%% done\n" else printf "\n" endif flush list REG = const X if opts.tvar REG += dummify(time) endif ols y REG --quiet BICi = BICrit ($ess, N, T, 1, K, s2hat) BIC |= {1, BICi, $ess} if BICi < bestBIC bestg = 1 endif printf "%s\n", "-"*32 printf " G\t BIC\t SSR\n" printf "%s\n", "-"*32 loop i = 1 .. G printf "%2d\t%7.4g\t%7.4g", i, BIC[G-i+1,2], BIC[G-i+1,3] if i == bestg printf " *" endif printf "\n" endloop printf "%s\n", "-"*32 opts.BICselect = 1 opts.verbose = verbose else mod = gfe_estimate_fixed(y, X, G, opts) endif err = add_sandwich(y, X, &mod) mod.etime = $stopwatch if opts.verbose printf "estimation done (obj = %g)\n", mod.obj endif mod.depvarname = argname(y) return mod pn = mod.parnames cs = mod.coeff ~ mod.stderr opts = mod.opts tvstr = opts.tvar ? "time-varying" : "time-invariant" printf "\n%s\n", "="*64 printf "Group fixed effects estimator " printf "(algorithm %d)\n", opts.algorithm if opts.algorithm == 1 printf "method = " if opts.method == "random" printf "%d random starts\n", opts.randstarts else printf "%s\n", opts.method endif else printf "%d random starts, ", opts.randstarts string extra = opts.parallel == 0 ? "single-threaded" : sprintf("%d processes", opts.parallel) printf "%s\n", extra endif printf "dependent variable: %s\n", mod.depvarname printf "Using %d units and %d time periods (%d total observations)\n", mod.N, mod.T, mod.N*mod.T printf "%d %s groups ", opts.ngrp, tvstr if mod.bootrep > 0 printf "(bootstrap s.e., %d replications)\n", mod.bootrep else printf "(sandwich s.e.)\n" endif printf "%s\n", "="*64 modprint cs pn printf "SSR = %g, Log-likelihood = %g\n", mod.obj, mod.lnl printf "Elapsed time = %g seconds\n\n", mod.etime if showfreq printf "Number of units per group:\n" showfreq(mod.grp) printf "\n" endif matrix paths = mshape(NA, mod.T, mod.N) G = mod.opts.ngrp if mod.opts.tvar loop i = 1 .. G sel = which(mod.grp .= i) paths[,sel] = mshape(mod.a[,i], $pd, nelem(sel)) endloop series ret = vec(paths) else tmp = replace(mod.grp, seq(1, G), mod.a) ret = pexpand(tmp) endif discrete ret setinfo ret --description="Group effect" return ret opts = mod.opts G = opts.ngrp out = exists(dest) ? dest : "display" tvar = opts.tvar matrix plt = mod.a s = seq(1, G)' lab = "group " ~ s title = "Grouped fixed effects" if !tvar # a is a column vector plt = plt ~ s rnameset(plt, lab) cnameset(plt, "alpha group") plot plt option with-boxes literal set linetype 1 pt 1 lc rgb "#888888" literal set boxwidth 0.75 relative literal set style fill solid 0.5 literal set title '@title' printf "set xrange [0:%d]", G+1 end plot --output=@out else plt = plt ~ seq(1, rows(plt))' lab += "time" cnameset(plt, lab) plot plt options with-lp single-yaxis literal set xlabel "Time" literal set ylabel "Group effect" literal set title '@title' printf "set xrange [0:%d]", $pd+1 literal set key outside right end plot --output=@out endif bundle ret = empty ret.ngrp = G ret.tvar = 1 ret.method = "random" ret.randstarts = NA ret.bootiter = 128 ret.algorithm = 2 ret.verbose = 1 ret.BICselect = 0 # use MPI? si = $sysinfo ret.parallel = si.mpi * si.ncores return ret list ALL = y || X series o = ok(ALL) return min(o) == 1 NT = N*T ret = ssr/NT + s2Gmax * (G*T + N + K)/NT * log(NT) return ret /* -doc- provides a grouping of the units into n groups, according to the parameter "method" key in the bundle b: random: discrete uniform given: basically a no-op, it just uses the matrix initgrp residuals: kmeans on pmeaned OLS residuals of y on X kmeans: kmeans on pmeaned y and X kmeans-res: kmeans on the NxT matrix of OLS residuals sorting: based on sorted OLS residuals of y on X The resulting grouping is stored into the bundle b, under the key "inigrp" */ startmth = b.method n = b.ngrp if startmth == "random" grp = mrandgen(i, 1, n, N, 1) elif startmth == "given" grp = b.inigrp elif startmth == "residuals" ols y const X --quiet mean_u = pshrink(pmean($uhat)) kmod = kmeans(mean_u, n) grp = kmod.cluster_id elif startmth == "kmeans-res" matrix U = {} mols({y}, 1 ~ {X}, &U) matrix U = mshape(U, $pd, N)' kmod = kmeans(U, n) grp = kmod.cluster_id elif startmth == "kmeans" matrix tmp = pshrink(pmean(y)) loop foreach i X tmp ~= pshrink(pmean(X.$i)) endloop kmod = kmeans(tmp, n) grp = kmod.cluster_id elif startmth == "sorting" ols y const X --quiet mean_u = pmean($uhat) m = msortby(pshrink(mean_u), 1) h = -msortby( -(m ~ diff(m)), 2 ) limits = h[1:n, 1] grp = 1 + sumr(m .> limits') else funcerr "unknown starting method\n" endif b.inigrp = grp /* -doc- Bonhomme-Manresa (2015) -- Updating step: compute the coefficients given the groupings */ # for the time-varying model, interaction of unit and time \alpha_{g,t} series group_id = timevar ? (pexpand(grp-1)) * P2 + time : pexpand(grp) scalar nreg = nelem(X) matrix e = {} matrix a = {} b = fols({y}, {X}, {group_id}, &e, , &a) obj = e'e loglik = -$nobs/2 * (2.8378770664 + log(obj/$nobs)) bundle ret = _(b, loglik, obj) ret.a = timevar ? mshape(a, $pd, ngrp) : a chk = rows(ret.b) + nelem(ret.a) - nreg ret.drop = chk < ngrp return ret /* -doc - compute the ssr per unit for all possible groupings */ ssr_mat = zeros(N, ngrp) if !timevar ng = rows(a) # should be equal to ngrp, but it may be not loop i = 1 .. ng tmp = pshrink(psum((res-a[i])^2)) ssr_mat[,i] = tmp endloop else loop i = 1 .. ngrp series tmp = mshape(a[,i], $nobs, 1) e = res - tmp ssr_mat[,i] = pshrink(psum(e^2)) endloop endif return ssr_mat /* -doc- Bonhomme-Manresa (2015) -- Assignment step: given the coefficients, assign the units to the various groups (also check that no group remains empty */ ssr_mat = ssr_by_unit(res, a, N, ngrp, timevar) # for each unit, find the group that minimises the ssr matrix new = iminr(ssr_mat) # sanity check if rows(values(new)) == 1 printf "What???? ng = %d\n", ng endif return new /* -doc- This function performs algorithm 1 with different group initialisation: "given" expects an initial input */ N = $nobs/$pd # number of units bundle ret = _(N=N, parnames=varnames(X_in)) X = X_in - const nreg = nelem(X) bundle ret = ret + _(initgrp=initgrp, drop=0) done = 0 iter = 0 # P2 is the smallest power of 2 strictly larger that $pd. For # example, if $pd = 3, then P2 = 4; if $pd = 4, P2 = 8, etc. # it is used with time-varying GFE to build interactions # (unused with time-fixed effects) scalar P2 = timevar ? 2^ceil(log2($pd+1)) : NA grp = initgrp loop while !done && iter < maxiter iter++ bundle tmpbun = update_step(y, X, ngrp, timevar, grp, P2) if tmpbun.drop ret.drop = 1 break else b = tmpbun.b a = tmpbun.a loglik = tmpbun.loglik obj = tmpbun.obj endif series res = y - lincomb(X, b) new = assignment_step(res, a, N, ngrp, timevar) # check if any individual has moved from one group to another eq = (grp .= new)' grp = new done = minr(eq) if verbose printf "iter = %d, loglik = %g, ", iter, loglik printf "switchers = %d", N - sum(eq) printf "\n" flush endif endloop if ret.drop == 0 ret.coeff = b ret.a = a ret.nreg = nreg # ret.iter = iter ret.grp = grp ret.lnl = loglik ret.obj = obj endif if verbose printf "Time: %g (%3d iterations), llik = %12.5f\n", $stopwatch, iter, loglik endif return ret /* -doc- performs estimation via algorithm 1 */ X = X_in - const T = $pd N = $nobs/T timevar = opts.tvar ngrp = opts.ngrp method = opts.method verbose = opts.verbose ret = _(opts=opts, N=N, T=T, X=X_in, depvarname=argname(y)) if method == "random" maxit = opts.randstarts bestcrit = $huge bundle best i = 0 dropped = 0 loop while i < maxit inigrp = mrandgen(i, 1, ngrp, N, 1) b = algo1_core(y, X, ngrp, timevar, inigrp,, verbose > 1) if b.drop == 1 dropped++ else i++ if b.obj < bestcrit bestcrit = b.obj best = b nbest = i if verbose printf "\rmodel %4d: %15.8f\n", i, bestcrit flush endif endif endif endloop ret = ret + best ret.dropped = dropped else b = algo1_core(y, X, ngrp, timevar, opts.inigrp,, verbose > 1) ret = ret + b endif # series ret.depvar = y string ret.depvarname = argname(y) list ret.X = X return ret scalar err = 0 matrix grp = mod.grp scalar timevar = mod.opts.tvar scalar P2 = timevar ? 2^ceil(log2($pd+1)) : NA series id = $unit scalar k = nelem(X) series group_id = timevar ? (pexpand(grp-1)) * P2 + time : pexpand(grp) ols y const X dummify(group_id) --cluster=id --quiet mod.vcv = $vcv[2:k+1, 2:k+1] mod.stderr = sqrt(diag(mod.vcv)) mod.bootrep = 0 return err scalar N = $nobs/$pd opts.ngrp = G init_groups(y, X, N, &opts) scalar verbose = opts.verbose if opts.algorithm == 2 bundle GOAT = _(obj=$huge) scalar nmax = NA scalar timevar = opts.tvar scalar P2 = timevar ? 2^ceil(log2($pd+1)) : NA randstart = opts.randstarts if verbose string extra = opts.parallel == 0 ? "single-thread" : "MPI-parallel" printf "Starting algorithm 2... (%s)\n", extra printf "random starts = %d\n", randstart flush endif if opts.parallel == 0 # single-threaded estimation loop outside_iter = 1 .. randstart best = algo2_estim(y, X, G, timevar, opts.inigrp, nmax) if verbose printf "Random start %d, obj = %.6g (best so far = %.6g)\n", outside_iter, best.obj, GOAT.obj flush endif if best.obj < GOAT.obj GOAT=best endif endloop opts.inigrp = GOAT.gamma else # rather arbitrary auto initialisation of realloc_max # if not user-given best = algo2_mpiblock(y, X, G, timevar, opts.parallel, nmax, randstart, verbose) opts.inigrp = best.gamma endif opts.method="given" # temporary endif mod = algo1_estim(y, X, opts) if opts.algorithm == 2 opts.method="random" endif mod.depvarname = argname(y) return mod /* -doc- needed for the bootstrap */ m = mshape({x}, T, N) series ret = vec(m[,resmat]) return ret series y = mod.depvar list X = mod.X err = 0 scalar k = nelem(X) matrix tgrp = mod.grp scalar timevar = mod.tvar scalar N = mod.N scalar ngroups = mod.ngrp matrix BB = mshape(NA, bootrep, k) # will hold the bootstrapped coefficients scalar iter = 0 scalar good = 0 scalar dropped = 0 scalar maxiter = 2^12 scalar maxit_for_boot = 128 matrix lliks = zeros(bootrep,1) bootopts = mod.opts bootopts.method = "random" bootopts.trials = maxit_for_boot bootopts.verbose = 0 bootopts.algorithm = mod.algorithm if bootopts.algorithm > 1 # common stuff for the two version of algorithm 2 bundle GOAT = _(obj=$huge) scalar nmax = NA timevar = bootopts.tvar G = bootopts.ngrp verbose = bootopts.verbose randstart = 8 if verbose printf "Starting algorithm 2...\n\n" flush endif endif printf "Starting bootstrap iterations\n" flush loop while (iter < bootrep) iter++ if iter%5 == 0 printf "boot %4d/%4d\r", iter, bootrep flush endif boot_draw = mrandgen(i, 1, N, 1, N) # resample the dep. var series by = panel_resample(y, N, $pd, boot_draw) # resample the explanatory variables list bX = empty loop foreach j X series b$j = panel_resample(X.$j, N, $pd, boot_draw) bX += b$j endloop # re-estimate init_groups(by, bX, N, &bootopts) if bootopts.algorithm == 2 loop outside_iter = 1 .. randstart best = algo2_estim(by, bX, G, timevar, bootopts.inigrp, nmax) if best.obj < GOAT.obj GOAT=best endif endloop bootopts.method="given" bootopts.inigrp = GOAT.gamma elif bootopts.algorithm == 3 # rather arbitrary auto initialisation of realloc_max # if not user-given scalar neigh_max = ok(nmax) ? nmax : floor(sqrt(N)) scalar j_max = 128 scalar P2 = timevar ? 2^ceil(log2($pd+1)) : NA best = algo2_mpiblock(by, bX, G, timevar, nmax, randstart) bootopts.method = "given" bootopts.inigrp = best.gamma endif bb = algo1_estim(y, X, bootopts) # printf "bb.loglik = %g\n", bb.lnl # printf "coeff =\n%12.6f\n", bb.coeff' lliks[iter] = bb.lnl BB[iter,] = bb.coeff' endloop vcv = mcov(BB) if 1 bias = meanc(BB)' - mod.coeff vcv += bias * bias' endif mod.bootrep = bootrep mod.vcv = vcv mod.stderr = sqrt(diag(vcv)) return err frq = sumr(values(grp) .= grp') rlab = "Group " ~ seq(1, rows(frq)) ~ ":" rnameset(frq, rlab) printf "\n%4d", frq flush /* -doc- here we take an intial group assignment (from 1 to ngrp) and produce a new assignment where n units are switched randomly to another group (modular artithmetic is used for optimal perfomance) */ matrix ret = grp scalar ngrp = max(grp) scalar N = rows(grp) # extract units without replacement matrix tmp = seq(1,N)' ~ muniform(N,1) matrix to_rep = msortby(tmp,2)[1:n,1] # draw new group matrix grep = grp[to_rep] matrix shuffle = mrandgen(i, 1, ngrp-1, n, 1) ret[to_rep] = ((grep + shuffle - 1) % ngrp) + 1 return ret /* -doc- Algorithm 2 by Bonhomme and Manresa (2015) */ scalar N = $nobs/$pd # number of units bundle ret = _(N=N, parnames=varnames(X_in), tvar=timevar, ngrp=G) X = X_in - const nreg = nelem(X) # mX = {X} # go with an initial grouping matrix grp = initgrp bundle ret = ret + _(initgrp=grp, drop=0) # rather arbitrary auto initialisation of neigh_max # if not user-given scalar neigh_max = ok(nmax) ? nmax : floor(sqrt(N)) scalar j_max = 128 scalar P2 = timevar ? 2^ceil(log2($pd+1)) : NA bundle best = _(obj=$huge) is_start_ok = 0 loop while !is_start_ok matrix grp = mrandgen(i, 1, G, N, 1) bundle tmp = update_step(y, X, G, timevar, grp, P2) series res = y - lincomb(X, tmp.b) matrix gprime = assignment_step(res, tmp.a, N, G, timevar) is_start_ok = rows(values(gprime)) == G endloop loop j = 0 .. j_max n = 1 loop while n < neigh_max matrix gg = relocate(gprime, n) bundle tmp = algo1_core(y, X, G, timevar, gg) if tmp.drop continue endif series res = y - lincomb(X, tmp.coeff) matrix newgrp = tmp.grp matrix candidate = ssr_by_unit(res, tmp.a, N, G, timevar) scalar newmin = sumc(minr(candidate)) if newmin < best.obj # printf "\tnewmin = %g (iter = %d)\n", newmin, j gprime = iminr(candidate) best.obj = newmin best.gamma = gprime else n++ endif flush endloop endloop return best scalar h = floor(rands/n) ret = mshape(h, 1, n) top_up = rands - sum(ret) if top_up > 0 ret[1:top_up] += 1 endif return ret err = 0 bundle to_mpi = _(G, timevar, verbose, nmax) string to_mpi.inimeth = "random" matrix to_mpi.Xl = X matrix to_mpi.randstarts = randstarts scalar to_mpi.seed = $seed bwrite(to_mpi, "parpass.xml", 1) act_cores = ncores > randstarts ? randstarts : ncores mpi --send-functions --send-data --quiet --np=act_cores bundle from_outside = bread("parpass.xml", 1) scalar G = from_outside.G scalar timevar = from_outside.timevar scalar verbose = from_outside.verbose scalar nmax = from_outside.nmax scalar randstarts = from_outside.randstarts scalar seed = from_outside.seed string inimeth = from_outside.inimeth list X = from_outside.Xl my_seed = seed + 32*$mpirank set seed my_seed starts = starts_per_core(randstarts, $mpisize) local_starts = starts[$mpirank+1] bestcrit = $huge loop i = 1 .. local_starts bundle new = algo2_estim(y, X, G, timevar, nmax) if verbose > 1 printf "Process %d: start %d/%d, obj = %.6g\n", $mpirank, i, local_starts, new.obj flush endif if new.obj < bestcrit bestcrit = new.obj best = new endif endloop mpibarrier() if $mpirank > 0 mpisend(best, 0) else bundles mods = array($mpisize) mods[1] = best loop i = 2 .. $mpisize mods[i] = mpirecv(i-1) endloop endif bestobj = $huge if $mpirank == 0 if verbose printf "\n\n" endif nbest = 1 loop i = 1 .. $mpisize if verbose printf "thread %03d, obj = %12.6f\n", i, mods[i].obj endif if mods[i].obj < bestobj nbest = i bestobj = mods[i].obj endif endloop bundle ret = mods[nbest] bwrite(ret, "best.xml", 1) endif end mpi GOAT=bread("best.xml", 1) return GOAT set verbose off include GFE.gfn open abdata.gdt --quiet set seed 8273647 list X = w k ys ### make panel balanced list ALL = n X smpl (YEAR > 1976) --restrict series d = ok(ALL) smpl pmin(d) == 1 --restrict ### do GFE estimation nGrp = 6 opts = _(verbose = 1) ### uncomment below to enable BIC-based model selection # opts.BICselect = 1 mod = gfe_estimate(n, X, nGrp, opts) gfe_printout(mod) ### use algorithm 1 (should give a worse solution in terms of SSR) opts.algorithm = 1 opts.randstarts = 2048 mod = gfe_estimate(n, X, nGrp, opts) gfe_printout(mod)