//
//  Session.swift
//  Metrika
//
//  Author:  Simon-Pierre Boucher
//  Contact: contact@spboucher.ai
//  Copyright © 2026 Simon-Pierre Boucher. All rights reserved.
//

import Foundation
import ZQData
import ZQGPU
import ZQGraphics
import ZQParser
import ZQPlanner
import ZQPlugins
import ZQStats

public struct ZQEngineError: Error, Equatable, Sendable, CustomStringConvertible {
    public var message: String
    public init(_ message: String) { self.message = message }
    public var description: String { message }
}

/// One interactive session: working dataset, RNG state, panel settings,
/// logging. Every command mutation — console, do-file, or plugin — goes
/// through `execute(_:)` so all frontends share one execution path
/// (CLAUDE.md §7).
public actor ZQSession {
    public private(set) var frame: ZQDataFrame
    public private(set) var lastPlot: ZQPlotSpec?
    public private(set) var seed: UInt64 = 123_456_789
    public private(set) var panelVariable: String?
    public private(set) var timeVariable: String?

    private let store: ZQDataStore
    private let parser: ZQCommandParser
    private let planner = ZQPlanner(gpuAvailable: ZQGPUBootstrap.isAvailable)
    private let scriptCommands: [String: ZQScriptCommand]
    private let pluginRegistry: ZQPluginRegistry
    private var logFileURL: URL?
    private var scriptDepth = 0
    private var lastEstimation: EstimationState?

    /// - Parameters:
    ///   - discoverUserCommands: scan the commands directory for `.zyq`
    ///     script commands.
    ///   - commandsDirectory: override of the discovery location (tests).
    ///   - plugins: SPM-compiled-in native plugins to register. A plugin
    ///     whose verb shadows a built-in command is rejected.
    public init(
        discoverUserCommands: Bool = true,
        commandsDirectory: URL? = nil,
        plugins: [any ZQCommandPlugin] = []
    ) throws {
        self.store = try ZQDataStore()
        self.frame = try ZQDataFrame()

        var commands: [String: ZQScriptCommand] = [:]
        if discoverUserCommands {
            let scripts = commandsDirectory.map {
                ZQPluginDiscovery.scriptCommands(in: $0)
            } ?? ZQPluginDiscovery.scriptCommands()
            for script in scripts {
                commands[script.verb] = script
            }
        }
        self.scriptCommands = commands

        let builtinVerbs = Set(ZQVerbTable.builtin.verbs.keys)
        self.pluginRegistry = try ZQPluginRegistry(
            plugins: plugins, reservedVerbs: builtinVerbs
        )

        // Plugin verbs join the grammar so their commands parse with full
        // varlist/if/in/options structure (no abbreviations for plugins).
        var verbs = ZQVerbTable.builtin.verbs
        for verb in pluginRegistry.verbs {
            verbs[verb] = verb.count
        }
        self.parser = ZQCommandParser(verbTable: ZQVerbTable(
            verbs: verbs,
            fileVerbs: ZQVerbTable.builtin.fileVerbs,
            assignmentVerbs: ZQVerbTable.builtin.assignmentVerbs,
            prefixVerbs: ZQVerbTable.builtin.prefixVerbs,
            compoundVerbs: ZQVerbTable.builtin.compoundVerbs
        ))
    }

    // MARK: - Entry point

    public func execute(_ line: String) async throws -> ZQResult {
        let result: ZQResult
        do {
            result = try await run(line)
        } catch {
            try? appendToLog(command: line, output: "error: \(error)")
            throw error
        }
        try? appendToLog(command: line, output: result.text)
        return result
    }

    /// Executes a do-file: runs each line, concatenating output. Stops at
    /// the first error, reporting the failing line number.
    public func executeScript(_ text: String) async throws -> ZQResult {
        var outputs: [String] = []
        var scalars: [String: Double] = [:]
        for (index, rawLine) in text.split(
            separator: "\n", omittingEmptySubsequences: false
        ).enumerated() {
            let line = String(rawLine)
            do {
                let result = try await execute(line)
                if !result.text.isEmpty {
                    outputs.append(". \(line)\n\(result.text)")
                }
                scalars.merge(result.scalars) { _, new in new }
            } catch {
                throw ZQEngineError("line \(index + 1): \(error)")
            }
        }
        return ZQResult(text: outputs.joined(separator: "\n\n"), scalars: scalars)
    }

    private func run(_ line: String) async throws -> ZQResult {
        // Script commands are routed BEFORE the parser: their verbs are
        // not in the grammar and their arguments are raw macro text.
        let trimmed = line.trimmingCharacters(in: .whitespaces)
        let verbWord = String(trimmed.prefix { !$0.isWhitespace })
        if let script = scriptCommands[verbWord] {
            guard scriptDepth < 8 else {
                throw ZQEngineError("user command recursion too deep")
            }
            scriptDepth += 1
            defer { scriptDepth -= 1 }
            let arguments = String(trimmed.dropFirst(verbWord.count))
                .trimmingCharacters(in: .whitespaces)
            let body = script.expandedLines(arguments: arguments)
                .joined(separator: "\n")
            return try await executeScript(body)
        }

        guard let command = try parser.parse(line) else {
            return ZQResult(text: "")
        }

        if let registered = pluginRegistry[command.verb] {
            return try await executePlugin(registered, command: command)
        }

        let plan = planner.plan(command, rowCount: frame.rowCount)
        return try await dispatch(plan.command, backend: plan.backend)
    }

    /// Runs a native plugin: validates the parsed command against its
    /// declared syntax, hands it a value copy of the dataset, and installs
    /// a returned replacement frame only when `mutates: true`.
    private func executePlugin(
        _ registered: ZQPluginRegistry.Registered, command: ZQCommand
    ) async throws -> ZQResult {
        if let violation = pluginRegistry.validate(command, against: registered.syntax) {
            throw ZQEngineError("\(registered.verb): \(violation)")
        }
        let context = ZQContext(command: command, frame: frame)
        var result = try await registered.plugin.execute(context)
        if let replacement = result.replacementFrame {
            guard registered.syntax.mutates else {
                throw ZQEngineError(
                    "\(registered.verb): plugin returned a dataset but does not declare 'mutates'"
                )
            }
            frame = replacement
            result.replacementFrame = nil
        }
        return result
    }

    private func dispatch(
        _ command: ZQCommand, backend: ZQBackend = .cpu
    ) async throws -> ZQResult {
        switch command.verb {
        case "use": return try await handleUse(command)
        case "sysuse": return try await handleSysuse(command)
        case "save": return try await handleSave(command)
        case "clear": return handleClear()
        case "describe": return handleDescribe()
        case "summarize": return try handleSummarize(command)
        case "generate": return try handleGenerate(command, replace: false)
        case "replace": return try handleGenerate(command, replace: true)
        case "drop": return try handleDropKeep(command, keep: false)
        case "keep": return try handleDropKeep(command, keep: true)
        case "count": return try handleCount(command)
        case "list": return try handleList(command)
        case "regress": return try handleRegress(command)
        case "xtreg": return try handleXTReg(command)
        case "ivregress": return try handleIVRegress(command)
        case "predict": return try handlePredict(command)
        case "margins": return try handleMargins(command)
        case "elasticnet": return try handleElasticNet(command, defaultAlpha: nil)
        case "boost": return try handleBoost(command)
        case "lasso": return try handleElasticNet(command, defaultAlpha: 1)
        case "logit": return try handleGLM(command, family: .logit)
        case "probit": return try handleGLM(command, family: .probit)
        case "poisson": return try handleGLM(command, family: .poisson)
        case "tabulate": return try handleTabulate(command)
        case "correlate": return try handleCorrelate(command)
        case "set": return try handleSet(command)
        case "xtset": return try handleXTSet(command)
        case "display": return try handleDisplay(command)
        case "bootstrap": return try await handleBootstrap(command, backend: backend)
        case "permute": return try await handlePermute(command)
        case "bayes": return try handleBayes(command)
        case "graph": return try handleGraph(command)
        case "histogram", "scatter", "kdensity":
            var promoted = command
            promoted.subverb = command.verb
            promoted.verb = "graph"
            return try handleGraph(promoted)
        case "log": return try handleLog(command)
        case "help": return handleHelp(command)
        default:
            throw ZQEngineError("command '\(command.verb)' is not implemented yet")
        }
    }

    // MARK: - Data commands

    private func handleUse(_ command: ZQCommand) async throws -> ZQResult {
        guard let path = command.argument else {
            throw ZQEngineError("use: file name required")
        }
        let url = URL(fileURLWithPath: (path as NSString).expandingTildeInPath)
        frame = try await store.load(contentsOf: url)
        return ZQResult(
            text: "(\(frame.rowCount) observations, \(frame.columns.count) variables loaded from \(url.lastPathComponent))",
            scalars: ["N": Double(frame.rowCount), "k": Double(frame.columns.count)]
        )
    }

    /// `sysuse name` — load a bundled sample dataset; bare `sysuse` lists
    /// what ships with the app.
    private func handleSysuse(_ command: ZQCommand) async throws -> ZQResult {
        let samples = ["sales", "mtcars"]
        guard let name = command.argument else {
            return ZQResult(
                text: "Sample datasets: \(samples.joined(separator: ", ")) — load one with 'sysuse <name>'"
            )
        }
        let stem = name.hasSuffix(".csv") ? String(name.dropLast(4)) : name
        guard samples.contains(stem),
              let url = Bundle.module.url(
                forResource: stem, withExtension: "csv", subdirectory: "Samples"
              ) else {
            throw ZQEngineError(
                "sysuse: unknown sample '\(name)' — available: \(samples.joined(separator: ", "))"
            )
        }
        frame = try await store.load(contentsOf: url)
        lastEstimation = nil
        return ZQResult(
            text: "(\(frame.rowCount) observations, \(frame.columns.count) variables — sample dataset '\(stem)')",
            scalars: ["N": Double(frame.rowCount), "k": Double(frame.columns.count)]
        )
    }

    private func handleSave(_ command: ZQCommand) async throws -> ZQResult {
        guard let path = command.argument else {
            throw ZQEngineError("save: file name required")
        }
        guard !frame.isEmpty else { throw ZQEngineError("no data in memory") }
        let url = URL(fileURLWithPath: (path as NSString).expandingTildeInPath)
        if FileManager.default.fileExists(atPath: url.path),
           !command.hasOption("replace") {
            throw ZQEngineError("file \(url.lastPathComponent) already exists — add ', replace' to overwrite")
        }
        try await store.save(frame, to: url)
        return ZQResult(text: "file \(url.lastPathComponent) saved")
    }

    private func handleClear() -> ZQResult {
        frame = try! ZQDataFrame()
        panelVariable = nil
        timeVariable = nil
        return ZQResult(text: "")
    }

    private func handleDescribe() -> ZQResult {
        guard !frame.isEmpty else {
            return ZQResult(text: "Contains data: 0 observations")
        }
        var lines = [
            "Contains data: \(frame.rowCount) observations, \(frame.columns.count) variables",
            "",
            TableFormatter.pad("Variable", 12) + "   " +
            TableFormatter.pad("Type", 8, right: false) + "  Missing",
            String(repeating: "-", count: 34),
        ]
        for column in frame.columns {
            let type = column.data.isNumeric ? "float64" : "string"
            lines.append(
                TableFormatter.pad(column.name, 12) + "   " +
                TableFormatter.pad(type, 8, right: false) + "  " +
                TableFormatter.pad("\(column.data.missingCount)", 7)
            )
        }
        return ZQResult(
            text: lines.joined(separator: "\n"),
            scalars: ["N": Double(frame.rowCount), "k": Double(frame.columns.count)]
        )
    }

    private func handleSummarize(_ command: ZQCommand) throws -> ZQResult {
        guard !frame.isEmpty else { throw ZQEngineError("no data in memory") }
        let mask = try observationMask(command)
        let detail = command.hasOption("detail")

        var names: [String] = command.varlist.flatMap(\.referencedNames)
        if names.isEmpty {
            names = frame.columns.filter(\.data.isNumeric).map(\.name)
        }

        var lines = [
            TableFormatter.pad("Variable", 12) + " | " +
            TableFormatter.pad("Obs", 9) + "  " +
            TableFormatter.pad("Mean", 11) + "  " +
            TableFormatter.pad("Std. dev.", 11) + "  " +
            TableFormatter.pad("Min", 10) + "  " +
            TableFormatter.pad("Max", 10),
            String(repeating: "-", count: 13) + "+" + String(repeating: "-", count: 62),
        ]
        var scalars: [String: Double] = [:]

        for name in names {
            let (values, missing) = try frame.requireNumeric(name)
            let restrictedMissing = (0..<values.count).map { missing[$0] || !mask[$0] }
            let summary = ZQSummarize.summary(
                name: name, values: values, missing: restrictedMissing, detail: detail
            )
            if let summary {
                lines.append(
                    TableFormatter.pad(name, 12) + " | " +
                    TableFormatter.pad("\(summary.observationCount)", 9) + "  " +
                    TableFormatter.pad(TableFormatter.general(summary.mean), 11) + "  " +
                    TableFormatter.pad(TableFormatter.general(summary.standardDeviation), 11) + "  " +
                    TableFormatter.pad(TableFormatter.general(summary.minimum), 10) + "  " +
                    TableFormatter.pad(TableFormatter.general(summary.maximum), 10)
                )
                if let detailBlock = summary.detail {
                    for p in [1, 5, 10, 25, 50, 75, 90, 95, 99] {
                        lines.append(
                            TableFormatter.pad("p\(p)", 16) + "  " +
                            TableFormatter.pad(
                                TableFormatter.general(detailBlock.percentiles[p] ?? .nan), 11
                            )
                        )
                    }
                    lines.append(
                        TableFormatter.pad("skewness", 16) + "  " +
                        TableFormatter.pad(TableFormatter.general(detailBlock.skewness), 11)
                    )
                    lines.append(
                        TableFormatter.pad("kurtosis", 16) + "  " +
                        TableFormatter.pad(TableFormatter.general(detailBlock.kurtosis), 11)
                    )
                }
                scalars["N"] = Double(summary.observationCount)
                scalars["mean"] = summary.mean
                scalars["sd"] = summary.standardDeviation
                scalars["min"] = summary.minimum
                scalars["max"] = summary.maximum
                scalars["Var"] = summary.variance
            } else {
                lines.append(
                    TableFormatter.pad(name, 12) + " | " +
                    TableFormatter.pad("0", 9)
                )
            }
        }
        return ZQResult(text: lines.joined(separator: "\n"), scalars: scalars)
    }

    private func handleGenerate(_ command: ZQCommand, replace: Bool) throws -> ZQResult {
        guard let assignment = command.assignment else {
            throw ZQEngineError("\(replace ? "replace" : "generate"): syntax is 'newvar = expression'")
        }
        guard !frame.isEmpty else { throw ZQEngineError("no data in memory") }

        let evaluator = ExpressionEvaluator(frame: frame)
        var (values, missing) = try evaluator.evaluateNumericColumn(assignment.expression)
        let mask = try observationMask(command)

        if replace {
            let (oldValues, oldMissing) = try frame.requireNumeric(assignment.target)
            var changed = 0
            for i in 0..<frame.rowCount {
                if mask[i] {
                    let differs = oldMissing[i] != missing[i]
                        || (!missing[i] && oldValues[i] != values[i])
                    if differs { changed += 1 }
                } else {
                    values[i] = oldValues[i]
                    missing[i] = oldMissing[i]
                }
            }
            try frame.replaceColumn(ZQColumn(
                name: assignment.target,
                data: .float64(values: values, missing: missing)
            ))
            return ZQResult(text: "(\(changed) real changes made)")
        } else {
            var missingCount = 0
            for i in 0..<frame.rowCount where !mask[i] {
                values[i] = .nan
                missing[i] = true
            }
            for flag in missing where flag { missingCount += 1 }
            try frame.addColumn(ZQColumn(
                name: assignment.target,
                data: .float64(values: values, missing: missing)
            ))
            let note = missingCount > 0 ? "(\(missingCount) missing values generated)" : ""
            return ZQResult(text: note)
        }
    }

    private func handleDropKeep(_ command: ZQCommand, keep: Bool) throws -> ZQResult {
        guard !frame.isEmpty else { throw ZQEngineError("no data in memory") }

        // Observation form: `drop if expr` / `keep if expr` / `keep in 1/10`.
        if command.varlist.isEmpty, command.condition != nil || command.range != nil {
            let mask = try observationMask(command)
            let keepMask = keep ? mask : mask.map { !$0 }
            let before = frame.rowCount
            frame = try frame.filtered(by: keepMask)
            return ZQResult(text: "(\(before - frame.rowCount) observations deleted)")
        }

        // Variable form: `drop x y` / `keep x y`.
        let names = command.varlist.flatMap(\.referencedNames)
        guard !names.isEmpty else {
            throw ZQEngineError("\(keep ? "keep" : "drop"): specify variables or an 'if' condition")
        }
        if keep {
            let dropped = frame.columnNames.filter { !names.contains($0) }
            try frame.dropColumns(dropped)
        } else {
            try frame.dropColumns(names)
        }
        return ZQResult(text: "")
    }

    private func handleCount(_ command: ZQCommand) throws -> ZQResult {
        guard !frame.isEmpty else { throw ZQEngineError("no data in memory") }
        let mask = try observationMask(command)
        let count = mask.count { $0 }
        return ZQResult(text: "  \(count)", scalars: ["N": Double(count)])
    }

    private func handleList(_ command: ZQCommand) throws -> ZQResult {
        guard !frame.isEmpty else { throw ZQEngineError("no data in memory") }
        let mask = try observationMask(command)
        var names = command.varlist.flatMap(\.referencedNames)
        if names.isEmpty { names = frame.columnNames }
        let columns = try names.map { try frame.requireColumn($0) }

        let rows = (0..<frame.rowCount).filter { mask[$0] }.prefix(50)
        var lines = [names.map { TableFormatter.pad($0, 12) }.joined(separator: "  ")]
        for row in rows {
            let cells = columns.map { column -> String in
                switch column.data {
                case .float64(let values, let missing):
                    return TableFormatter.pad(
                        missing[row] ? "." : TableFormatter.general(values[row]), 12
                    )
                case .string(let values):
                    return TableFormatter.pad(values[row] ?? "", 12)
                }
            }
            lines.append(cells.joined(separator: "  "))
        }
        return ZQResult(text: lines.joined(separator: "\n"))
    }

    // MARK: - Settings

    private func handleSet(_ command: ZQCommand) throws -> ZQResult {
        switch command.subverb {
        case "seed":
            guard let text = command.argument, let value = UInt64(text) else {
                throw ZQEngineError("set seed: integer seed required")
            }
            seed = value
            return ZQResult(text: "")
        default:
            throw ZQEngineError("set: unknown setting '\(command.subverb ?? "")'")
        }
    }

    private func handleXTSet(_ command: ZQCommand) throws -> ZQResult {
        let names = command.varlist.flatMap(\.referencedNames)
        guard (1...2).contains(names.count) else {
            throw ZQEngineError("xtset: syntax is 'xtset panelvar [timevar]'")
        }
        for name in names {
            _ = try frame.requireColumn(name)
        }
        panelVariable = names[0]
        timeVariable = names.count > 1 ? names[1] : nil
        var text = "Panel variable: \(names[0])"
        if let timeVariable { text += "\nTime variable: \(timeVariable)" }
        return ZQResult(text: text)
    }

    private func handleDisplay(_ command: ZQCommand) throws -> ZQResult {
        guard let assignment = command.assignment else {
            throw ZQEngineError("display: expression required")
        }
        let evaluator = ExpressionEvaluator(frame: frame)
        switch try evaluator.evaluate(assignment.expression) {
        case .numericScalar(let value):
            return ZQResult(text: TableFormatter.general(value, significant: 10))
        case .stringScalar(let value):
            return ZQResult(text: value)
        case .missingScalar:
            return ZQResult(text: ".")
        case .numeric(let values, let missing):
            guard let first = values.first else { return ZQResult(text: ".") }
            return ZQResult(
                text: missing[0] ? "." : TableFormatter.general(first, significant: 10)
            )
        case .strings(let values):
            return ZQResult(text: values.first.flatMap { $0 } ?? "")
        }
    }

    private func handleLog(_ command: ZQCommand) throws -> ZQResult {
        if command.argument == "close" || command.subverb == "close" {
            logFileURL = nil
            return ZQResult(text: "(log closed)")
        }
        guard var path = command.argument else {
            throw ZQEngineError("log: syntax is 'log using filename' or 'log close'")
        }
        if path == "close" {
            logFileURL = nil
            return ZQResult(text: "(log closed)")
        }
        if path.hasPrefix("using ") { path = String(path.dropFirst(6)) }
        let url = URL(fileURLWithPath: (path as NSString).expandingTildeInPath)
        if !FileManager.default.fileExists(atPath: url.path) {
            FileManager.default.createFile(atPath: url.path, contents: nil)
        }
        logFileURL = url
        return ZQResult(text: "(log started: \(url.path))")
    }

    /// `help [command]` — renders the shared command reference.
    private func handleHelp(_ command: ZQCommand) -> ZQResult {
        let requested = command.varlist.flatMap(\.referencedNames).first
            ?? command.argument

        if let requested {
            guard let doc = ZQCommandReference.doc(for: requested) else {
                return ZQResult(text: "help: no entry for '\(requested)' — type 'help' for the full list")
            }
            var lines = [
                "\(doc.verb)\(doc.abbreviation.map { " (abbreviation: \($0))" } ?? "")",
                String(repeating: "-", count: 60),
                doc.summary,
                "",
                "Syntax:  \(doc.syntax)",
            ]
            if !doc.options.isEmpty {
                lines.append("")
                lines.append("Options:")
                for option in doc.options {
                    lines.append("  " + TableFormatter.pad(option.name, 16, right: false) + option.meaning)
                }
            }
            if !doc.examples.isEmpty {
                lines.append("")
                lines.append("Examples:")
                for example in doc.examples { lines.append("  . \(example)") }
            }
            if let notes = doc.notes {
                lines.append("")
                lines.append(notes)
            }
            return ZQResult(text: lines.joined(separator: "\n"))
        }

        var lines = ["Metrika command reference — 'help <command>' for details", ""]
        for category in ZQCommandReference.categories {
            let docs = ZQCommandReference.all.filter { $0.category == category }
            guard !docs.isEmpty else { continue }
            lines.append(category)
            for doc in docs {
                lines.append(
                    "  " + TableFormatter.pad(doc.verb, 14, right: false) + doc.summary.prefix(60)
                )
            }
            lines.append("")
        }
        return ZQResult(text: lines.joined(separator: "\n"))
    }

    private func appendToLog(command: String, output: String) throws {
        guard let logFileURL else { return }
        let entry = ". \(command)\n\(output)\n\n"
        let handle = try FileHandle(forWritingTo: logFileURL)
        defer { try? handle.close() }
        try handle.seekToEnd()
        try handle.write(contentsOf: Data(entry.utf8))
    }

    // MARK: - Estimation

    /// How a fitted regressor is recomputed from the working dataset —
    /// the estimation-state piece that lets `predict` (and later
    /// `margins`) evaluate xb on arbitrary observations.
    enum PredictorDefinition: Equatable, Sendable {
        case column(String)
        case indicator(variable: String, level: Double)
        case product([String])
        case constant
    }

    /// Stata's e() analog: what the last estimation command fitted.
    struct EstimationState: Sendable {
        enum Kind: Equatable, Sendable {
            case ols, iv
            case glm(ZQGLMFamily)
            case boost(ZQBoostModel)
        }

        var kind: Kind
        var responseName: String
        /// Coefficients aligned with their recomputation recipes.
        var coefficients: [(name: String, value: Double, definition: PredictorDefinition)]
        /// Covariance matrix, column-major k×k aligned with `coefficients`.
        var vce: [Double]
        /// Degrees of freedom for t inference; nil means normal (ML).
        var inferenceDF: Double?
        /// Estimation-sample design (column-major n×k incl. constant) —
        /// kept for GLMs, whose average marginal effects need it.
        var sampleDesign: [Double]?
        var sampleSize: Int
    }

    /// Assembled estimation sample after if/in restriction and listwise
    /// deletion (§6: explicit report line for dropped observations).
    struct RegressionSample {
        var responseName: String
        var y: [Double]
        var predictors: [(name: String, values: [Double])]
        /// Aligned with `predictors`.
        var definitions: [PredictorDefinition]
        var clusterLabels: [Int]?
        var droppedMissing: Int
    }

    private func buildRegressionSample(
        _ command: ZQCommand, clusterVariable: String?
    ) throws -> RegressionSample {
        guard !frame.isEmpty else { throw ZQEngineError("no data in memory") }
        guard let first = command.varlist.first else {
            throw ZQEngineError("\(command.verb): dependent variable required")
        }
        guard case .simple(let responseName) = first else {
            throw ZQEngineError("\(command.verb): dependent variable cannot use factor notation")
        }

        let mask = try observationMask(command)
        let (yAll, yMissing) = try frame.requireNumeric(responseName)

        // Collect raw regressor sources (numeric columns) for missingness.
        struct Source {
            var spec: ZQVarSpec
            var columns: [(name: String, values: [Double], missing: [Bool])]
        }
        var sources: [Source] = []
        for spec in command.varlist.dropFirst() {
            var columns: [(String, [Double], [Bool])] = []
            for name in spec.referencedNames {
                let (values, missing) = try frame.requireNumeric(name)
                columns.append((name, values, missing))
            }
            sources.append(Source(spec: spec, columns: columns))
        }

        var clusterColumn: ZQColumn?
        if let clusterVariable {
            clusterColumn = try frame.requireColumn(clusterVariable)
        }

        // Listwise deletion over response, regressors, and cluster var.
        var keptRows: [Int] = []
        var droppedMissing = 0
        for i in 0..<frame.rowCount where mask[i] {
            var complete = !yMissing[i]
            if complete {
                outer: for source in sources {
                    for column in source.columns where column.missing[i] {
                        complete = false
                        break outer
                    }
                }
            }
            if complete, let clusterColumn {
                switch clusterColumn.data {
                case .float64(_, let missing): complete = !missing[i]
                case .string(let values): complete = values[i] != nil
                }
            }
            if complete {
                keptRows.append(i)
            } else {
                droppedMissing += 1
            }
        }
        guard !keptRows.isEmpty else {
            throw ZQEngineError("no observations")
        }

        let y = keptRows.map { yAll[$0] }

        // Expand specs on the estimation sample (plan-time expansion),
        // recording each column's recomputation recipe for predict.
        var predictors: [(name: String, values: [Double])] = []
        var definitions: [PredictorDefinition] = []
        for source in sources {
            switch source.spec {
            case .simple(let name):
                let (values, _) = try frame.requireNumeric(name)
                predictors.append((name, keptRows.map { values[$0] }))
                definitions.append(.column(name))
            case .factor(let op, let name):
                let (values, missing) = try frame.requireNumeric(name)
                switch op {
                case "i", "ib":
                    let sampleValues = keptRows.map { values[$0] }
                    let sampleMissing = keptRows.map { missing[$0] }
                    let expanded = ZQFactorExpansion.expand(
                        name: name, values: sampleValues, missing: sampleMissing
                    )
                    guard !expanded.isEmpty else {
                        throw ZQEngineError("factor variable '\(name)' has a single level")
                    }
                    for term in expanded {
                        predictors.append((term.name, term.values))
                        definitions.append(.indicator(variable: name, level: term.level))
                    }
                case "c":
                    predictors.append((name, keptRows.map { values[$0] }))
                    definitions.append(.column(name))
                default:
                    throw ZQEngineError("unsupported factor operator '\(op).'")
                }
            case .interaction(let parts, _):
                // Continuous-by-continuous interactions only.
                var label: [String] = []
                var factors: [String] = []
                var product = [Double](repeating: 1, count: keptRows.count)
                for part in parts {
                    guard case .factor(let op, let name) = part, op == "c" else {
                        guard case .simple(let name) = part else {
                            throw ZQEngineError("factor interactions beyond c.#c. are not supported yet")
                        }
                        let (values, _) = try frame.requireNumeric(name)
                        for (j, row) in keptRows.enumerated() { product[j] *= values[row] }
                        label.append(name)
                        factors.append(name)
                        continue
                    }
                    let (values, _) = try frame.requireNumeric(name)
                    for (j, row) in keptRows.enumerated() { product[j] *= values[row] }
                    label.append("c.\(name)")
                    factors.append(name)
                }
                predictors.append((label.joined(separator: "#"), product))
                definitions.append(.product(factors))
            }
        }

        // Cluster labels → dense integer codes.
        var clusterLabels: [Int]?
        if let clusterColumn {
            var codes: [Int] = []
            codes.reserveCapacity(keptRows.count)
            switch clusterColumn.data {
            case .float64(let values, _):
                var mapping: [Double: Int] = [:]
                for row in keptRows {
                    let value = values[row]
                    if let code = mapping[value] {
                        codes.append(code)
                    } else {
                        mapping[value] = mapping.count
                        codes.append(mapping.count - 1)
                    }
                }
            case .string(let values):
                var mapping: [String: Int] = [:]
                for row in keptRows {
                    let value = values[row] ?? ""
                    if let code = mapping[value] {
                        codes.append(code)
                    } else {
                        mapping[value] = mapping.count
                        codes.append(mapping.count - 1)
                    }
                }
            }
            clusterLabels = codes
        }

        return RegressionSample(
            responseName: responseName,
            y: y,
            predictors: predictors,
            definitions: definitions,
            clusterLabels: clusterLabels,
            droppedMissing: droppedMissing
        )
    }

    /// Records the estimation state (Stata's e()) for post-estimation
    /// commands. The constant, when fitted, is the last coefficient.
    private func recordEstimation(
        kind: EstimationState.Kind,
        responseName: String,
        coefficients: [ZQCoefficient],
        definitions: [PredictorDefinition],
        includeConstant: Bool,
        vce: [Double] = [],
        inferenceDF: Double? = nil,
        sampleDesign: [Double]? = nil,
        sampleSize: Int = 0
    ) {
        var recipes = definitions
        if includeConstant { recipes.append(.constant) }
        guard recipes.count == coefficients.count else {
            lastEstimation = nil
            return
        }
        lastEstimation = EstimationState(
            kind: kind,
            responseName: responseName,
            coefficients: zip(coefficients, recipes).map {
                (name: $0.name, value: $0.estimate, definition: $1)
            },
            vce: vce,
            inferenceDF: inferenceDF,
            sampleDesign: sampleDesign,
            sampleSize: sampleSize
        )
    }

    /// Column-major estimation-sample design from sample predictors.
    private func designMatrix(
        from sample: RegressionSample, includeConstant: Bool
    ) -> [Double] {
        let n = sample.y.count
        var design = [Double]()
        design.reserveCapacity(n * (sample.predictors.count + 1))
        for column in sample.predictors { design.append(contentsOf: column.values) }
        if includeConstant { design.append(contentsOf: [Double](repeating: 1, count: n)) }
        return design
    }

    private func varianceEstimator(
        _ command: ZQCommand, clusterLabels: [Int]?
    ) throws -> ZQVarianceEstimator {
        if let clusterLabels { return .cluster(clusterLabels) }
        if command.hasOption("hc2") { return .hc2 }
        if command.hasOption("hc3") { return .hc3 }
        if command.hasOption("robust") || command.hasOption("r") { return .hc1 }
        if let vce = command.option("vce")?.firstArgument {
            switch vce {
            case "robust": return .hc1
            case "hc2": return .hc2
            case "hc3": return .hc3
            default:
                throw ZQEngineError("vce(\(vce)) is not supported")
            }
        }
        return .classical
    }

    private func handleRegress(_ command: ZQCommand) throws -> ZQResult {
        let clusterVariable = command.option("cluster")?.firstArgument
            ?? clusterFromVCE(command)
        let sample = try buildRegressionSample(command, clusterVariable: clusterVariable)
        let variance = try varianceEstimator(command, clusterLabels: sample.clusterLabels)
        let level = try confidenceLevel(command)

        let result = try ZQOLS.fit(
            y: sample.y,
            predictors: sample.predictors,
            includeConstant: !command.hasOption("noconstant"),
            variance: variance,
            confidenceLevel: level
        )
        recordEstimation(
            kind: .ols,
            responseName: sample.responseName,
            coefficients: result.coefficients,
            definitions: sample.definitions,
            includeConstant: !command.hasOption("noconstant"),
            vce: result.vce,
            inferenceDF: result.inferenceDF,
            sampleSize: result.observationCount
        )

        var header: [String] = []
        switch variance {
        case .classical: header.append("Linear regression")
        case .cluster: header.append("Linear regression (cluster-robust)")
        default: header.append("Linear regression (robust)")
        }
        if sample.droppedMissing > 0 {
            header.append("(\(sample.droppedMissing) observations dropped due to missing values)")
        }

        var stats = [
            ("Number of obs", TableFormatter.general(Double(result.observationCount))),
            ("R-squared", TableFormatter.fixed(result.rSquared, decimals: 4)),
            ("Adj R-squared", TableFormatter.fixed(result.adjustedRSquared, decimals: 4)),
            ("Root MSE", TableFormatter.general(result.rootMSE, significant: 6)),
        ]
        if let f = result.fStatistic, let df = result.fDF, let p = result.fPValue {
            stats.insert(
                ("F(\(Int(df.0)), \(Int(df.1)))", TableFormatter.fixed(f, decimals: 2)),
                at: 1
            )
            stats.insert(("Prob > F", TableFormatter.fixed(p, decimals: 4)), at: 2)
        }
        if let g = result.clusterCount {
            stats.append(("Clusters", "\(g)"))
        }
        for (label, value) in stats {
            header.append(
                TableFormatter.pad(label, 46, right: false) + "= " +
                TableFormatter.pad(value, 12)
            )
        }

        var lines = header
        lines.append("")
        let widths = [12, 12, 11, 8, 8, 22]
        lines.append(
            TableFormatter.pad(sample.responseName, widths[0]) + " | " +
            TableFormatter.pad("Coefficient", widths[1]) + "  " +
            TableFormatter.pad("Std. err.", widths[2]) + "  " +
            TableFormatter.pad("t", widths[3]) + "  " +
            TableFormatter.pad("P>|t|", widths[4]) + "  " +
            TableFormatter.pad("[\(Int(level * 100))% conf. interval]", widths[5])
        )
        lines.append(TableFormatter.rule(widths))
        var scalars: [String: Double] = [
            "N": Double(result.observationCount),
            "r2": result.rSquared,
            "r2_a": result.adjustedRSquared,
            "rmse": result.rootMSE,
            "df_r": result.inferenceDF,
        ]
        if let f = result.fStatistic { scalars["F"] = f }
        if let g = result.clusterCount { scalars["N_clust"] = Double(g) }

        for coefficient in result.coefficients {
            lines.append(
                TableFormatter.pad(coefficient.name, widths[0]) + " | " +
                TableFormatter.pad(TableFormatter.general(coefficient.estimate), widths[1]) + "  " +
                TableFormatter.pad(TableFormatter.general(coefficient.standardError), widths[2]) + "  " +
                TableFormatter.pad(TableFormatter.fixed(coefficient.tStatistic, decimals: 2), widths[3]) + "  " +
                TableFormatter.pad(TableFormatter.fixed(coefficient.pValue, decimals: 3), widths[4]) + "  " +
                TableFormatter.pad(
                    TableFormatter.general(coefficient.confidenceLower) + "   " +
                    TableFormatter.general(coefficient.confidenceUpper),
                    widths[5]
                )
            )
            scalars["b_\(coefficient.name)"] = coefficient.estimate
            scalars["se_\(coefficient.name)"] = coefficient.standardError
        }
        return ZQResult(text: lines.joined(separator: "\n"), scalars: scalars)
    }

    private func handleXTReg(_ command: ZQCommand) throws -> ZQResult {
        guard command.hasOption("fe") else {
            throw ZQEngineError("xtreg: only the fixed-effects estimator is implemented — add ', fe'")
        }
        guard let panelVariable else {
            throw ZQEngineError("xtreg: declare the panel first with 'xtset panelvar [timevar]'")
        }
        // v0.2 restriction: the cluster variable must be the panel variable.
        var clustered = false
        if let clusterName = command.option("cluster")?.firstArgument
            ?? clusterFromVCE(command) {
            guard clusterName == panelVariable else {
                throw ZQEngineError(
                    "xtreg: clustering on a variable other than the panel (\(panelVariable)) is not supported yet"
                )
            }
            clustered = true
        }

        // Panel codes ride the same listwise-deletion path as cluster codes.
        let sample = try buildRegressionSample(command, clusterVariable: panelVariable)
        guard let panelCodes = sample.clusterLabels else {
            throw ZQEngineError("xtreg: panel variable missing from sample")
        }

        let result = try ZQFixedEffects.fitWithin(
            y: sample.y,
            predictors: sample.predictors,
            groups: panelCodes,
            clustered: clustered,
            confidenceLevel: try confidenceLevel(command)
        )
        // predict after xtreg needs the estimated u_i — not stored yet.
        lastEstimation = nil

        var header = ["Fixed-effects (within) regression"]
        if sample.droppedMissing > 0 {
            header.append("(\(sample.droppedMissing) observations dropped due to missing values)")
        }
        var stats: [(String, String)] = [
            ("Number of obs", "\(result.observationCount)"),
            ("Number of groups", "\(result.groupCount)"),
            ("R-squared (within)", TableFormatter.fixed(result.rSquaredWithin, decimals: 4)),
            ("Root MSE", TableFormatter.general(result.rootMSE, significant: 6)),
        ]
        if let f = result.fStatistic, let p = result.fPValue {
            stats.append(("F", TableFormatter.fixed(f, decimals: 2)))
            stats.append(("Prob > F", TableFormatter.fixed(p, decimals: 4)))
        }
        for (label, value) in stats {
            header.append(
                TableFormatter.pad(label, 46, right: false) + "= " +
                TableFormatter.pad(value, 12)
            )
        }

        var scalars: [String: Double] = [
            "N": Double(result.observationCount),
            "N_g": Double(result.groupCount),
            "r2_w": result.rSquaredWithin,
            "rmse": result.rootMSE,
            "df_r": result.inferenceDF,
        ]
        let lines = header + [""] + coefficientTable(
            response: sample.responseName,
            coefficients: result.coefficients,
            statisticLabel: "t",
            level: try confidenceLevel(command),
            scalars: &scalars
        )
        return ZQResult(text: lines.joined(separator: "\n"), scalars: scalars)
    }

    private func handleIVRegress(_ command: ZQCommand) throws -> ZQResult {
        guard command.subverb == "2sls" else {
            throw ZQEngineError("ivregress: only the 2sls estimator is implemented")
        }
        guard let ivSpec = command.ivSpec else {
            throw ZQEngineError("ivregress: syntax is 'ivregress 2sls depvar [exog] (endog = instruments)'")
        }
        guard !frame.isEmpty else { throw ZQEngineError("no data in memory") }
        guard let first = command.varlist.first, case .simple(let response) = first else {
            throw ZQEngineError("ivregress: dependent variable required")
        }
        let exogenousNames = command.varlist.dropFirst().flatMap(\.referencedNames)
        let clusterVariable = command.option("cluster")?.firstArgument
            ?? clusterFromVCE(command)

        // Listwise deletion across every variable involved.
        let allNames = [response] + exogenousNames + ivSpec.endogenous + ivSpec.instruments
        let sources = try allNames.map { try frame.requireNumeric($0) }
        var clusterColumn: ZQColumn?
        if let clusterVariable {
            clusterColumn = try frame.requireColumn(clusterVariable)
        }
        let mask = try observationMask(command)
        var keptRows: [Int] = []
        var dropped = 0
        for i in 0..<frame.rowCount where mask[i] {
            var complete = sources.allSatisfy { !$0.missing[i] }
            if complete, let clusterColumn {
                switch clusterColumn.data {
                case .float64(_, let missing): complete = !missing[i]
                case .string(let values): complete = values[i] != nil
                }
            }
            if complete { keptRows.append(i) } else { dropped += 1 }
        }
        guard !keptRows.isEmpty else { throw ZQEngineError("no observations") }

        func gather(_ names: [String]) throws -> [(name: String, values: [Double])] {
            try names.map { name in
                let (values, _) = try frame.requireNumeric(name)
                return (name, keptRows.map { values[$0] })
            }
        }
        let y = try frame.requireNumeric(response).values
        let clusterLabels: [Int]? = try clusterVariable.map { name in
            let column = try frame.requireColumn(name)
            var mapping: [String: Int] = [:]
            var codes: [Int] = []
            for row in keptRows {
                let key: String
                switch column.data {
                case .float64(let values, _): key = "\(values[row])"
                case .string(let values): key = values[row] ?? ""
                }
                if let code = mapping[key] {
                    codes.append(code)
                } else {
                    mapping[key] = mapping.count
                    codes.append(mapping.count - 1)
                }
            }
            return codes
        }

        let result = try ZQIV.fit2SLS(
            y: keptRows.map { y[$0] },
            endogenous: try gather(ivSpec.endogenous),
            exogenous: try gather(exogenousNames),
            instruments: try gather(ivSpec.instruments),
            includeConstant: !command.hasOption("noconstant"),
            variance: try varianceEstimator(command, clusterLabels: clusterLabels),
            confidenceLevel: try confidenceLevel(command)
        )
        recordEstimation(
            kind: .iv,
            responseName: response,
            coefficients: result.coefficients,
            definitions: (ivSpec.endogenous + exogenousNames).map { .column($0) },
            includeConstant: !command.hasOption("noconstant"),
            vce: result.vce,
            inferenceDF: result.inferenceDF,
            sampleSize: result.observationCount
        )

        var header = ["Instrumental variables (2SLS) regression"]
        if dropped > 0 {
            header.append("(\(dropped) observations dropped due to missing values)")
        }
        header.append("Instrumented: \(ivSpec.endogenous.joined(separator: " "))")
        header.append("Instruments:  \(ivSpec.instruments.joined(separator: " "))")
        var stats: [(String, String)] = [
            ("Number of obs", "\(result.observationCount)"),
            ("R-squared", TableFormatter.fixed(result.rSquared, decimals: 4)),
            ("Root MSE", TableFormatter.general(result.rootMSE, significant: 6)),
        ]
        if let f = result.fStatistic, let p = result.fPValue {
            stats.insert(("F", TableFormatter.fixed(f, decimals: 2)), at: 1)
            stats.insert(("Prob > F", TableFormatter.fixed(p, decimals: 4)), at: 2)
        }
        if let g = result.clusterCount { stats.append(("Clusters", "\(g)")) }
        for (label, value) in stats {
            header.append(
                TableFormatter.pad(label, 46, right: false) + "= " +
                TableFormatter.pad(value, 12)
            )
        }

        var scalars: [String: Double] = [
            "N": Double(result.observationCount),
            "r2": result.rSquared,
            "rmse": result.rootMSE,
            "df_r": result.inferenceDF,
        ]
        if let f = result.fStatistic { scalars["F"] = f }
        let lines = header + [""] + coefficientTable(
            response: response,
            coefficients: result.coefficients,
            statisticLabel: "t",
            level: try confidenceLevel(command),
            scalars: &scalars
        )
        return ZQResult(text: lines.joined(separator: "\n"), scalars: scalars)
    }

    /// Shared coefficient-table rendering (regress/xtreg/ivregress/GLMs).
    private func coefficientTable(
        response: String,
        coefficients: [ZQCoefficient],
        statisticLabel: String,
        level: Double,
        scalars: inout [String: Double]
    ) -> [String] {
        let widths = [12, 12, 11, 8, 8, 22]
        var lines = [
            TableFormatter.pad(response, widths[0]) + " | " +
            TableFormatter.pad("Coefficient", widths[1]) + "  " +
            TableFormatter.pad("Std. err.", widths[2]) + "  " +
            TableFormatter.pad(statisticLabel, widths[3]) + "  " +
            TableFormatter.pad("P>|\(statisticLabel)|", widths[4]) + "  " +
            TableFormatter.pad("[\(Int(level * 100))% conf. interval]", widths[5]),
            TableFormatter.rule(widths),
        ]
        for coefficient in coefficients {
            lines.append(
                TableFormatter.pad(coefficient.name, widths[0]) + " | " +
                TableFormatter.pad(TableFormatter.general(coefficient.estimate), widths[1]) + "  " +
                TableFormatter.pad(TableFormatter.general(coefficient.standardError), widths[2]) + "  " +
                TableFormatter.pad(TableFormatter.fixed(coefficient.tStatistic, decimals: 2), widths[3]) + "  " +
                TableFormatter.pad(TableFormatter.fixed(coefficient.pValue, decimals: 3), widths[4]) + "  " +
                TableFormatter.pad(
                    TableFormatter.general(coefficient.confidenceLower) + "   " +
                    TableFormatter.general(coefficient.confidenceUpper),
                    widths[5]
                )
            )
            scalars["b_\(coefficient.name)"] = coefficient.estimate
            scalars["se_\(coefficient.name)"] = coefficient.standardError
        }
        return lines
    }

    private func handleGLM(
        _ command: ZQCommand, family: ZQGLMFamily
    ) throws -> ZQResult {
        let clusterVariable = command.option("cluster")?.firstArgument
            ?? clusterFromVCE(command)
        let sample = try buildRegressionSample(command, clusterVariable: clusterVariable)
        let variance = try varianceEstimator(command, clusterLabels: sample.clusterLabels)
        let level = try confidenceLevel(command)

        let result = try ZQGLM.fit(
            y: sample.y,
            predictors: sample.predictors,
            family: family,
            includeConstant: !command.hasOption("noconstant"),
            variance: variance,
            confidenceLevel: level
        )
        recordEstimation(
            kind: .glm(family),
            responseName: sample.responseName,
            coefficients: result.coefficients,
            definitions: sample.definitions,
            includeConstant: !command.hasOption("noconstant"),
            vce: result.vce,
            sampleDesign: designMatrix(
                from: sample,
                includeConstant: !command.hasOption("noconstant")
            ),
            sampleSize: result.observationCount
        )

        let title: String
        switch family {
        case .logit: title = "Logistic regression"
        case .probit: title = "Probit regression"
        case .poisson: title = "Poisson regression"
        }

        var header = [title]
        if sample.droppedMissing > 0 {
            header.append("(\(sample.droppedMissing) observations dropped due to missing values)")
        }
        var stats: [(String, String)] = [
            ("Number of obs", "\(result.observationCount)"),
        ]
        if let chi2 = result.chiSquared, let p = result.chiSquaredPValue {
            let label: String
            if case .classical = variance {
                label = "LR chi2(\(result.chiSquaredDF))"
            } else {
                label = "Wald chi2(\(result.chiSquaredDF))"
            }
            stats.append((label, TableFormatter.fixed(chi2, decimals: 2)))
            stats.append(("Prob > chi2", TableFormatter.fixed(p, decimals: 4)))
        }
        stats.append(("Log likelihood", TableFormatter.fixed(result.logLikelihood, decimals: 4)))
        stats.append(("Pseudo R2", TableFormatter.fixed(result.pseudoRSquared, decimals: 4)))
        if let g = result.clusterCount {
            stats.append(("Clusters", "\(g)"))
        }
        for (label, value) in stats {
            header.append(
                TableFormatter.pad(label, 46, right: false) + "= " +
                TableFormatter.pad(value, 12)
            )
        }

        var lines = header
        lines.append("")
        let widths = [12, 12, 11, 8, 8, 22]
        lines.append(
            TableFormatter.pad(sample.responseName, widths[0]) + " | " +
            TableFormatter.pad("Coefficient", widths[1]) + "  " +
            TableFormatter.pad("Std. err.", widths[2]) + "  " +
            TableFormatter.pad("z", widths[3]) + "  " +
            TableFormatter.pad("P>|z|", widths[4]) + "  " +
            TableFormatter.pad("[\(Int(level * 100))% conf. interval]", widths[5])
        )
        lines.append(TableFormatter.rule(widths))

        var scalars: [String: Double] = [
            "N": Double(result.observationCount),
            "ll": result.logLikelihood,
            "ll_0": result.nullLogLikelihood,
            "r2_p": result.pseudoRSquared,
        ]
        if let chi2 = result.chiSquared { scalars["chi2"] = chi2 }
        if let g = result.clusterCount { scalars["N_clust"] = Double(g) }

        for coefficient in result.coefficients {
            lines.append(
                TableFormatter.pad(coefficient.name, widths[0]) + " | " +
                TableFormatter.pad(TableFormatter.general(coefficient.estimate), widths[1]) + "  " +
                TableFormatter.pad(TableFormatter.general(coefficient.standardError), widths[2]) + "  " +
                TableFormatter.pad(TableFormatter.fixed(coefficient.tStatistic, decimals: 2), widths[3]) + "  " +
                TableFormatter.pad(TableFormatter.fixed(coefficient.pValue, decimals: 3), widths[4]) + "  " +
                TableFormatter.pad(
                    TableFormatter.general(coefficient.confidenceLower) + "   " +
                    TableFormatter.general(coefficient.confidenceUpper),
                    widths[5]
                )
            )
            scalars["b_\(coefficient.name)"] = coefficient.estimate
            scalars["se_\(coefficient.name)"] = coefficient.standardError
        }
        return ZQResult(text: lines.joined(separator: "\n"), scalars: scalars)
    }

    // MARK: - Tabulate & correlate

    private func handleTabulate(_ command: ZQCommand) throws -> ZQResult {
        guard !frame.isEmpty else { throw ZQEngineError("no data in memory") }
        let names = command.varlist.flatMap(\.referencedNames)
        guard names.count == 1 || names.count == 2 else {
            throw ZQEngineError("tabulate: syntax is 'tabulate varname [varname2]'")
        }
        let mask = try observationMask(command)

        /// Rendered level label per observation; nil = missing (excluded
        /// unless the `missing` option is given).
        func labels(_ name: String) throws -> [String?] {
            let column = try frame.requireColumn(name)
            switch column.data {
            case .float64(let values, let missing):
                return (0..<frame.rowCount).map { i in
                    missing[i] ? nil : TableFormatter.general(values[i])
                }
            case .string(let values):
                return values
            }
        }

        let includeMissing = command.hasOption("missing")
        let first = try labels(names[0])

        if names.count == 1 {
            var counts: [String: Int] = [:]
            var total = 0
            for i in 0..<frame.rowCount where mask[i] {
                guard let label = includeMissing ? (first[i] ?? ".") : first[i] else {
                    continue
                }
                counts[label, default: 0] += 1
                total += 1
            }
            guard total > 0 else { throw ZQEngineError("no observations") }

            var lines = [
                TableFormatter.pad(names[0], 14) + " | " +
                TableFormatter.pad("Freq.", 9) + "  " +
                TableFormatter.pad("Percent", 9) + "  " +
                TableFormatter.pad("Cum.", 9),
                String(repeating: "-", count: 15) + "+" + String(repeating: "-", count: 33),
            ]
            var cumulative = 0.0
            for label in counts.keys.sorted(by: numericAwareLess) {
                let count = counts[label]!
                let percent = 100 * Double(count) / Double(total)
                cumulative += percent
                lines.append(
                    TableFormatter.pad(label, 14) + " | " +
                    TableFormatter.pad("\(count)", 9) + "  " +
                    TableFormatter.pad(TableFormatter.fixed(percent, decimals: 2), 9) + "  " +
                    TableFormatter.pad(TableFormatter.fixed(cumulative, decimals: 2), 9)
                )
            }
            lines.append(String(repeating: "-", count: 15) + "+" + String(repeating: "-", count: 33))
            lines.append(
                TableFormatter.pad("Total", 14) + " | " +
                TableFormatter.pad("\(total)", 9) + "  " +
                TableFormatter.pad("100.00", 9)
            )
            return ZQResult(
                text: lines.joined(separator: "\n"),
                scalars: ["N": Double(total), "r": Double(counts.count)]
            )
        }

        // Two-way table.
        let second = try labels(names[1])
        var cells: [String: [String: Int]] = [:]
        var rowTotals: [String: Int] = [:]
        var columnTotals: [String: Int] = [:]
        var total = 0
        for i in 0..<frame.rowCount where mask[i] {
            let rowLabel = includeMissing ? (first[i] ?? ".") : first[i]
            let columnLabel = includeMissing ? (second[i] ?? ".") : second[i]
            guard let rowLabel, let columnLabel else { continue }
            cells[rowLabel, default: [:]][columnLabel, default: 0] += 1
            rowTotals[rowLabel, default: 0] += 1
            columnTotals[columnLabel, default: 0] += 1
            total += 1
        }
        guard total > 0 else { throw ZQEngineError("no observations") }

        let rows = rowTotals.keys.sorted(by: numericAwareLess)
        let columns = columnTotals.keys.sorted(by: numericAwareLess)
        let width = 9
        var lines = [
            TableFormatter.pad(names[0], 14) + " | " +
            columns.map { TableFormatter.pad($0, width) }.joined(separator: " ") +
            "  " + TableFormatter.pad("Total", width),
            String(repeating: "-", count: 15) + "+" +
            String(repeating: "-", count: (width + 1) * (columns.count + 1) + 2),
        ]
        for row in rows {
            let cellsText = columns.map {
                TableFormatter.pad("\(cells[row]?[$0] ?? 0)", width)
            }.joined(separator: " ")
            lines.append(
                TableFormatter.pad(row, 14) + " | " + cellsText + "  " +
                TableFormatter.pad("\(rowTotals[row] ?? 0)", width)
            )
        }
        lines.append(
            TableFormatter.pad("Total", 14) + " | " +
            columns.map { TableFormatter.pad("\(columnTotals[$0] ?? 0)", width) }
                .joined(separator: " ") +
            "  " + TableFormatter.pad("\(total)", width)
        )
        return ZQResult(
            text: lines.joined(separator: "\n"),
            scalars: [
                "N": Double(total),
                "r": Double(rows.count),
                "c": Double(columns.count),
            ]
        )
    }

    /// Sorts numeric labels numerically, everything else lexically.
    private func numericAwareLess(_ a: String, _ b: String) -> Bool {
        if let x = Double(a), let y = Double(b) { return x < y }
        return a < b
    }

    private func handleCorrelate(_ command: ZQCommand) throws -> ZQResult {
        guard !frame.isEmpty else { throw ZQEngineError("no data in memory") }
        var names = command.varlist.flatMap(\.referencedNames)
        if names.isEmpty {
            names = frame.columns.filter(\.data.isNumeric).map(\.name)
        }
        guard names.count >= 2 else {
            throw ZQEngineError("correlate: at least two numeric variables required")
        }
        let mask = try observationMask(command)
        let sources = try names.map { try frame.requireNumeric($0) }

        // Listwise deletion across the varlist (Stata correlate).
        var kept: [Int] = []
        for i in 0..<frame.rowCount where mask[i] {
            if sources.allSatisfy({ !$0.missing[i] }) { kept.append(i) }
        }
        guard kept.count > 1 else { throw ZQEngineError("no observations") }

        let columns = sources.map { source in kept.map { source.values[$0] } }
        let matrix = ZQCorrelate.matrix(columns: columns)

        var lines = [
            "(obs=\(kept.count))",
            "",
            TableFormatter.pad("", 12) + " | " +
            names.map { TableFormatter.pad($0, 9) }.joined(separator: " "),
            String(repeating: "-", count: 13) + "+" +
            String(repeating: "-", count: (9 + 1) * names.count + 1),
        ]
        for (i, name) in names.enumerated() {
            // Lower triangle only, Stata-style.
            let cells = (0...i).map {
                TableFormatter.pad(TableFormatter.fixed(matrix[i][$0], decimals: 4), 9)
            }.joined(separator: " ")
            lines.append(TableFormatter.pad(name, 12) + " | " + cells)
        }
        return ZQResult(
            text: lines.joined(separator: "\n"),
            scalars: ["N": Double(kept.count), "rho": matrix[0][1]]
        )
    }

    private func clusterFromVCE(_ command: ZQCommand) -> String? {
        guard let vce = command.option("vce"),
              vce.arguments.count == 2,
              vce.arguments[0] == "cluster" else { return nil }
        return vce.arguments[1]
    }

    private func confidenceLevel(_ command: ZQCommand) throws -> Double {
        guard let text = command.option("level")?.firstArgument else { return 0.95 }
        guard let value = Double(text), value > 0, value < 100 else {
            throw ZQEngineError("level() must be between 0 and 100")
        }
        return value / 100
    }

    // MARK: - Post-estimation

    /// `predict newvar [, xb | residuals | pr | n]` — evaluates the last
    /// estimation over ALL current observations (missing inputs yield
    /// missing predictions). Defaults: xb after regress/ivregress, pr
    /// after logit/probit, n (mean count) after poisson.
    private func handlePredict(_ command: ZQCommand) throws -> ZQResult {
        guard let estimation = lastEstimation else {
            throw ZQEngineError("predict: no estimation results — run a regression first")
        }
        guard !frame.isEmpty else { throw ZQEngineError("no data in memory") }
        let names = command.varlist.flatMap(\.referencedNames)
        guard names.count == 1, let target = names.first else {
            throw ZQEngineError("predict: syntax is 'predict newvar [, statistic]'")
        }

        // Boosted models predict through their trees, not a linear xb.
        if case .boost(let model) = estimation.kind {
            return try predictBoost(
                model: model, target: target, command: command,
                responseName: estimation.responseName
            )
        }

        // Statistic selection and validation against the model kind.
        enum Statistic { case xb, residuals, probability, mean }
        var statistic: Statistic
        switch estimation.kind {
        case .ols, .iv: statistic = .xb
        case .glm(.logit), .glm(.probit): statistic = .probability
        case .glm(.poisson): statistic = .mean
        case .boost: fatalError("handled above")
        }
        for option in command.options {
            switch option.name {
            case "xb": statistic = .xb
            case "residuals", "resid": statistic = .residuals
            case "pr":
                switch estimation.kind {
                case .glm(.logit), .glm(.probit): statistic = .probability
                default:
                    throw ZQEngineError("predict: 'pr' requires logit or probit results")
                }
            case "n":
                guard case .glm(.poisson) = estimation.kind else {
                    throw ZQEngineError("predict: 'n' requires poisson results")
                }
                statistic = .mean
            default:
                throw ZQEngineError("predict: unknown statistic '\(option.name)'")
            }
        }

        // Linear predictor over the full dataset.
        let n = frame.rowCount
        var xb = [Double](repeating: 0, count: n)
        var missing = [Bool](repeating: false, count: n)
        for coefficient in estimation.coefficients {
            switch coefficient.definition {
            case .constant:
                for i in 0..<n { xb[i] += coefficient.value }
            case .column(let name):
                let (values, columnMissing) = try frame.requireNumeric(name)
                for i in 0..<n {
                    if columnMissing[i] { missing[i] = true } else {
                        xb[i] += coefficient.value * values[i]
                    }
                }
            case .indicator(let variable, let level):
                let (values, columnMissing) = try frame.requireNumeric(variable)
                for i in 0..<n {
                    if columnMissing[i] { missing[i] = true }
                    else if values[i] == level { xb[i] += coefficient.value }
                }
            case .product(let variables):
                let columns = try variables.map { try frame.requireNumeric($0) }
                for i in 0..<n {
                    var product = 1.0
                    for column in columns {
                        if column.missing[i] { missing[i] = true; break }
                        product *= column.values[i]
                    }
                    if !missing[i] { xb[i] += coefficient.value * product }
                }
            }
        }

        // Transform.
        var result = [Double](repeating: .nan, count: n)
        var label: String
        switch statistic {
        case .xb:
            label = "linear prediction"
            for i in 0..<n where !missing[i] { result[i] = xb[i] }
        case .residuals:
            label = "residuals"
            let (yValues, yMissing) = try frame.requireNumeric(estimation.responseName)
            for i in 0..<n where !missing[i] && !yMissing[i] {
                switch estimation.kind {
                case .ols, .iv:
                    result[i] = yValues[i] - xb[i]
                case .glm(.logit):
                    result[i] = yValues[i] - 1 / (1 + Foundation.exp(-xb[i]))
                case .glm(.probit):
                    result[i] = yValues[i] - ZQDistributions.normalCDF(xb[i])
                case .glm(.poisson):
                    result[i] = yValues[i] - Foundation.exp(xb[i])
                case .boost:
                    break   // handled by predictBoost
                }
            }
            for i in 0..<n where yMissing[i] { missing[i] = true }
        case .probability:
            label = "predicted probability"
            for i in 0..<n where !missing[i] {
                if case .glm(.probit) = estimation.kind {
                    result[i] = ZQDistributions.normalCDF(xb[i])
                } else {
                    result[i] = 1 / (1 + Foundation.exp(-xb[i]))
                }
            }
        case .mean:
            label = "predicted mean count"
            for i in 0..<n where !missing[i] { result[i] = Foundation.exp(xb[i]) }
        }
        for i in 0..<n where missing[i] { result[i] = .nan }

        try frame.addColumn(ZQColumn(
            name: target,
            data: .float64(values: result, missing: missing)
        ))
        let missingCount = missing.count { $0 }
        var note = "(\(label) '\(target)' generated"
        if missingCount > 0 { note += ", \(missingCount) missing values" }
        note += ")"
        return ZQResult(text: note, scalars: ["N": Double(n - missingCount)])
    }

    /// `elasticnet y x…, lambda(#) [alpha(#)]` and `lasso y x…, lambda(#)`
    /// — penalized linear regression, glmnet conventions. Without
    /// lambda(), reports λ_max (the smallest all-zero penalty) as a hint.
    private func handleElasticNet(
        _ command: ZQCommand, defaultAlpha: Double?
    ) throws -> ZQResult {
        let sample = try buildRegressionSample(command, clusterVariable: nil)

        var alpha = defaultAlpha ?? 1
        if let text = command.option("alpha")?.firstArgument {
            guard defaultAlpha == nil else {
                throw ZQEngineError("lasso: alpha is fixed at 1 — use elasticnet for other mixes")
            }
            guard let value = Double(text) else {
                throw ZQEngineError("elasticnet: invalid alpha()")
            }
            alpha = value
        }
        guard let lambdaText = command.option("lambda")?.firstArgument,
              let lambda = Double(lambdaText) else {
            let hint = ZQElasticNet.lambdaMax(
                y: sample.y, predictors: sample.predictors, alpha: alpha
            )
            throw ZQEngineError(
                "\(command.verb): lambda(#) required — lambda_max for this data is \(TableFormatter.general(hint, significant: 6))"
            )
        }

        let result = try ZQElasticNet.fit(
            y: sample.y,
            predictors: sample.predictors,
            alpha: alpha,
            lambda: lambda
        )

        // predict works after penalized fits; margins (needing a VCE)
        // deliberately does not.
        recordEstimation(
            kind: .ols,
            responseName: sample.responseName,
            coefficients: result.coefficients.map {
                ZQCoefficient(
                    name: $0.name, estimate: $0.value, standardError: .nan,
                    tStatistic: .nan, pValue: .nan,
                    confidenceLower: .nan, confidenceUpper: .nan
                )
            } + [ZQCoefficient(
                name: "_cons", estimate: result.intercept, standardError: .nan,
                tStatistic: .nan, pValue: .nan,
                confidenceLower: .nan, confidenceUpper: .nan
            )],
            definitions: sample.definitions,
            includeConstant: true,
            vce: [],
            sampleSize: sample.y.count
        )

        var lines = [
            alpha == 1 ? "Lasso linear model" : "Elastic-net linear model (alpha = \(TableFormatter.general(alpha)))",
        ]
        if sample.droppedMissing > 0 {
            lines.append("(\(sample.droppedMissing) observations dropped due to missing values)")
        }
        for (label, value) in [
            ("Number of obs", "\(sample.y.count)"),
            ("lambda", TableFormatter.general(lambda, significant: 6)),
            ("Nonzero coefficients", "\(result.nonzeroCount) of \(result.coefficients.count)"),
            ("R-squared", TableFormatter.fixed(result.rSquared, decimals: 4)),
        ] {
            lines.append(
                TableFormatter.pad(label, 46, right: false) + "= " +
                TableFormatter.pad(value, 12)
            )
        }
        lines.append("")
        lines.append(
            TableFormatter.pad(sample.responseName, 12) + " | " +
            TableFormatter.pad("Coefficient", 12)
        )
        lines.append(String(repeating: "-", count: 13) + "+" + String(repeating: "-", count: 14))

        var scalars: [String: Double] = [
            "N": Double(sample.y.count),
            "lambda": lambda,
            "alpha": alpha,
            "k_nonzero": Double(result.nonzeroCount),
            "r2": result.rSquared,
        ]
        for coefficient in result.coefficients where coefficient.value != 0 {
            lines.append(
                TableFormatter.pad(coefficient.name, 12) + " | " +
                TableFormatter.pad(TableFormatter.general(coefficient.value), 12)
            )
            scalars["b_\(coefficient.name)"] = coefficient.value
        }
        // Zeroed coefficients still surface as scalars (as exact 0).
        for coefficient in result.coefficients where coefficient.value == 0 {
            scalars["b_\(coefficient.name)"] = 0
        }
        lines.append(
            TableFormatter.pad("_cons", 12) + " | " +
            TableFormatter.pad(TableFormatter.general(result.intercept), 12)
        )
        scalars["b__cons"] = result.intercept
        return ZQResult(text: lines.joined(separator: "\n"), scalars: scalars)
    }

    /// `boost y x…, rounds(#) [eta(#) maxdepth(#) lambda(#)]` —
    /// gradient-boosted regression trees (xgboost exact-greedy, squared
    /// loss). Deterministic: no subsampling.
    private func handleBoost(_ command: ZQCommand) throws -> ZQResult {
        let sample = try buildRegressionSample(command, clusterVariable: nil)
        func numberOption(_ name: String, default defaultValue: Double) throws -> Double {
            guard let text = command.option(name)?.firstArgument else { return defaultValue }
            guard let value = Double(text), value > 0 else {
                throw ZQEngineError("boost: invalid \(name)()")
            }
            return value
        }
        guard let roundsText = command.option("rounds")?.firstArgument,
              let rounds = Int(roundsText), rounds > 0 else {
            throw ZQEngineError("boost: rounds(#) required")
        }
        let model = try ZQGradientBoosting.fit(
            y: sample.y,
            features: sample.predictors,
            rounds: rounds,
            learningRate: try numberOption("eta", default: 0.3),
            maxDepth: Int(try numberOption("maxdepth", default: 6)),
            lambda: try numberOption("lambda", default: 1)
        )

        // Training fit.
        let n = sample.y.count
        var rss = 0.0
        var tss = 0.0
        let meanY = sample.y.reduce(0, +) / Double(n)
        for i in 0..<n {
            let vector = sample.predictors.map { $0.values[i] }
            let predicted = model.predict(vector)
            rss += (sample.y[i] - predicted) * (sample.y[i] - predicted)
            tss += (sample.y[i] - meanY) * (sample.y[i] - meanY)
        }
        let rSquared = tss > 0 ? 1 - rss / tss : .nan
        let rmse = (rss / Double(n)).squareRoot()

        recordEstimation(
            kind: .boost(model),
            responseName: sample.responseName,
            coefficients: sample.predictors.map {
                ZQCoefficient(
                    name: $0.name, estimate: .nan, standardError: .nan,
                    tStatistic: .nan, pValue: .nan,
                    confidenceLower: .nan, confidenceUpper: .nan
                )
            },
            definitions: sample.definitions,
            includeConstant: false,
            sampleSize: n
        )

        var lines = ["Gradient-boosted trees (squared loss)"]
        if sample.droppedMissing > 0 {
            lines.append("(\(sample.droppedMissing) observations dropped due to missing values)")
        }
        for (label, value) in [
            ("Number of obs", "\(n)"),
            ("Trees", "\(rounds)"),
            ("Learning rate (eta)", TableFormatter.general(model.learningRate)),
            ("Max depth", command.option("maxdepth")?.firstArgument ?? "6"),
            ("Training R-squared", TableFormatter.fixed(rSquared, decimals: 4)),
            ("Training RMSE", TableFormatter.general(rmse, significant: 6)),
        ] {
            lines.append(
                TableFormatter.pad(label, 46, right: false) + "= " +
                TableFormatter.pad(value, 12)
            )
        }
        lines.append("")
        lines.append("Use 'predict newvar' for fitted values. Training fit is in-sample —")
        lines.append("expect it to be optimistic relative to held-out data.")
        return ZQResult(
            text: lines.joined(separator: "\n"),
            scalars: [
                "N": Double(n), "rounds": Double(rounds),
                "r2": rSquared, "rmse": rmse,
            ]
        )
    }

    /// `margins, dydx(varlist)` — average marginal effects with
    /// delta-method standard errors. OLS/IV effects are the coefficients
    /// themselves; GLM effects average dμ/dx over the estimation sample:
    ///   AME_j = β_j · (1/N) Σᵢ g′(xbᵢ)
    ///   ∂AME_j/∂β_m = (1/N) Σᵢ g″(xbᵢ)·x_im·β_j + δ_jm·(1/N) Σᵢ g′(xbᵢ)
    /// Only continuous (plain-column) regressors are supported so far.
    private func handleMargins(_ command: ZQCommand) throws -> ZQResult {
        guard let estimation = lastEstimation else {
            throw ZQEngineError("margins: no estimation results — run a regression first")
        }
        guard let dydx = command.option("dydx"), !dydx.arguments.isEmpty else {
            throw ZQEngineError("margins: syntax is 'margins, dydx(varlist)'")
        }
        guard !estimation.vce.isEmpty else {
            throw ZQEngineError("margins: the last estimation stored no covariance matrix")
        }

        let k = estimation.coefficients.count
        struct Effect {
            var name: String
            var dydx: Double
            var standardError: Double
        }
        var effects: [Effect] = []

        for name in dydx.arguments {
            guard let j = estimation.coefficients.firstIndex(where: {
                $0.definition == .column(name)
            }) else {
                if estimation.coefficients.contains(where: {
                    if case .indicator(let variable, _) = $0.definition { return variable == name }
                    if case .product(let variables) = $0.definition { return variables.contains(name) }
                    return false
                }) {
                    throw ZQEngineError(
                        "margins: dydx over factor or interaction terms is not supported yet"
                    )
                }
                throw ZQEngineError("margins: '\(name)' is not a regressor in the last estimation")
            }
            let beta = estimation.coefficients[j].value

            switch estimation.kind {
            case .boost:
                throw ZQEngineError("margins: not available after boost")
            case .ols, .iv:
                // Linear model: the marginal effect is the coefficient.
                effects.append(Effect(
                    name: name,
                    dydx: beta,
                    standardError: estimation.vce[j * k + j].squareRoot()
                ))

            case .glm(let family):
                guard let design = estimation.sampleDesign,
                      estimation.sampleSize > 0 else {
                    throw ZQEngineError("margins: estimation sample unavailable")
                }
                let n = estimation.sampleSize
                // xb over the estimation sample.
                var xb = [Double](repeating: 0, count: n)
                for m in 0..<k {
                    let value = estimation.coefficients[m].value
                    for i in 0..<n { xb[i] += value * design[m * n + i] }
                }
                // g′ and g″ per observation.
                var meanFirst = 0.0
                var second = [Double](repeating: 0, count: n)
                for i in 0..<n {
                    switch family {
                    case .logit:
                        let p = 1 / (1 + Foundation.exp(-xb[i]))
                        meanFirst += p * (1 - p)
                        second[i] = p * (1 - p) * (1 - 2 * p)
                    case .probit:
                        let density = Foundation.exp(-0.5 * xb[i] * xb[i])
                            / (2 * Double.pi).squareRoot()
                        meanFirst += density
                        second[i] = -xb[i] * density
                    case .poisson:
                        let mu = Foundation.exp(xb[i])
                        meanFirst += mu
                        second[i] = mu
                    }
                }
                meanFirst /= Double(n)

                // Delta-method gradient.
                var gradient = [Double](repeating: 0, count: k)
                for m in 0..<k {
                    var meanSecondX = 0.0
                    for i in 0..<n { meanSecondX += second[i] * design[m * n + i] }
                    gradient[m] = beta * meanSecondX / Double(n)
                }
                gradient[j] += meanFirst

                var variance = 0.0
                for a in 0..<k {
                    for b in 0..<k {
                        variance += gradient[a] * estimation.vce[b * k + a] * gradient[b]
                    }
                }
                effects.append(Effect(
                    name: name,
                    dydx: beta * meanFirst,
                    standardError: variance.squareRoot()
                ))
            }
        }

        // Render: t inference when the model carries a df, z otherwise.
        let level = try confidenceLevel(command)
        let usesT = estimation.inferenceDF != nil
        let critical: Double
        if let df = estimation.inferenceDF {
            critical = ZQDistributions.studentTQuantile(0.5 + level / 2, df: df)
        } else {
            critical = ZQDistributions.normalQuantile(0.5 + level / 2)
        }

        var lines = [
            "Average marginal effects                      Number of obs = " +
            TableFormatter.pad("\(estimation.sampleSize)", 10),
            "",
        ]
        let widths = [12, 12, 11, 8, 8, 22]
        let statLabel = usesT ? "t" : "z"
        lines.append(
            TableFormatter.pad("", widths[0]) + " | " +
            TableFormatter.pad("dy/dx", widths[1]) + "  " +
            TableFormatter.pad("Std. err.", widths[2]) + "  " +
            TableFormatter.pad(statLabel, widths[3]) + "  " +
            TableFormatter.pad("P>|\(statLabel)|", widths[4]) + "  " +
            TableFormatter.pad("[\(Int(level * 100))% conf. interval]", widths[5])
        )
        lines.append(TableFormatter.rule(widths))

        var scalars: [String: Double] = ["N": Double(estimation.sampleSize)]
        for effect in effects {
            let statistic = effect.dydx / effect.standardError
            let p: Double
            if let df = estimation.inferenceDF {
                p = ZQDistributions.tTestPValue(statistic, df: df)
            } else {
                p = 2 * (1 - ZQDistributions.normalCDF(abs(statistic)))
            }
            lines.append(
                TableFormatter.pad(effect.name, widths[0]) + " | " +
                TableFormatter.pad(TableFormatter.general(effect.dydx), widths[1]) + "  " +
                TableFormatter.pad(TableFormatter.general(effect.standardError), widths[2]) + "  " +
                TableFormatter.pad(TableFormatter.fixed(statistic, decimals: 2), widths[3]) + "  " +
                TableFormatter.pad(TableFormatter.fixed(p, decimals: 3), widths[4]) + "  " +
                TableFormatter.pad(
                    TableFormatter.general(effect.dydx - critical * effect.standardError) + "   " +
                    TableFormatter.general(effect.dydx + critical * effect.standardError),
                    widths[5]
                )
            )
            scalars["dydx_\(effect.name)"] = effect.dydx
            scalars["se_\(effect.name)"] = effect.standardError
        }
        return ZQResult(text: lines.joined(separator: "\n"), scalars: scalars)
    }

    /// Boosted-tree predictions: trees route missing features down the
    /// default (left) branch, so predictions exist for every observation.
    private func predictBoost(
        model: ZQBoostModel, target: String, command: ZQCommand,
        responseName: String
    ) throws -> ZQResult {
        for option in command.options
        where !["xb", "residuals", "resid"].contains(option.name) {
            throw ZQEngineError("predict: '\(option.name)' is not available after boost")
        }
        let wantsResiduals = command.hasOption("residuals") || command.hasOption("resid")

        let n = frame.rowCount
        let columns = try model.featureNames.map { try frame.requireNumeric($0) }
        var values = [Double](repeating: .nan, count: n)
        var missing = [Bool](repeating: false, count: n)
        var yValues = [Double](repeating: 0, count: n)
        var yMissing = [Bool](repeating: false, count: n)
        if wantsResiduals {
            (yValues, yMissing) = try frame.requireNumeric(responseName)
        }

        for i in 0..<n {
            let vector = columns.map { $0.missing[i] ? Double.nan : $0.values[i] }
            let predicted = model.predict(vector)
            if wantsResiduals {
                if yMissing[i] {
                    missing[i] = true
                } else {
                    values[i] = yValues[i] - predicted
                }
            } else {
                values[i] = predicted
            }
        }
        try frame.addColumn(ZQColumn(
            name: target, data: .float64(values: values, missing: missing)
        ))
        let label = wantsResiduals ? "residuals" : "boosted prediction"
        return ZQResult(text: "(\(label) '\(target)' generated)")
    }

    // MARK: - Bootstrap prefix

    private func handleBootstrap(
        _ command: ZQCommand, backend: ZQBackend
    ) async throws -> ZQResult {
        guard let body = command.body else {
            throw ZQEngineError("bootstrap: syntax is 'bootstrap, reps(#): command'")
        }
        guard body.verb == "regress" else {
            throw ZQEngineError("bootstrap currently supports 'regress' bodies only")
        }
        guard let repsText = command.option("reps")?.firstArgument,
              let reps = Int(repsText), reps >= 2 else {
            throw ZQEngineError("bootstrap: reps(#) with # ≥ 2 required")
        }
        if let seedText = command.option("seed")?.firstArgument {
            guard let value = UInt64(seedText) else {
                throw ZQEngineError("bootstrap: invalid seed")
            }
            seed = value
        }

        let sample = try buildRegressionSample(body, clusterVariable: nil)
        let includeConstant = !body.hasOption("noconstant")
        let point = try ZQOLS.fit(
            y: sample.y,
            predictors: sample.predictors,
            includeConstant: includeConstant,
            variance: .classical
        )

        let n = sample.y.count
        let k = point.coefficients.count
        let generator = Philox4x32(seed: seed)
        let y = sample.y
        let predictors = sample.predictors

        // Pairs bootstrap: replicate r regenerates its indices from the
        // Philox counter stream — deterministic and order-independent, so
        // replicates can run in parallel CPU chunks or as GPU batches with
        // identical resamples.
        let allDraws: [[Double]]
        if backend == .gpu {
            allDraws = ZQGPUBootstrap.pairsBootstrapOLS(
                y: y,
                predictors: predictors,
                includeConstant: includeConstant,
                replicates: reps,
                seed: seed
            )
        } else {
            allDraws = try await withThrowingTaskGroup(
                of: [(Int, [Double])].self
            ) { group in
                let chunkCount = min(
                    max(1, ProcessInfo.processInfo.activeProcessorCount), reps
                )
                let chunkSize = (reps + chunkCount - 1) / chunkCount
                for chunk in 0..<chunkCount {
                    let lower = chunk * chunkSize
                    let upper = min(reps, lower + chunkSize)
                    guard lower < upper else { continue }
                    group.addTask {
                        var results: [(Int, [Double])] = []
                        results.reserveCapacity(upper - lower)
                        for replicate in lower..<upper {
                            let indices = ZQResampling.pairsBootstrapIndices(
                                replicate: replicate, sampleSize: n, generator: generator
                            )
                            let resampledY = indices.map { y[$0] }
                            let resampledPredictors = predictors.map { column in
                                (name: column.name, values: indices.map { column.values[$0] })
                            }
                            let fit = try ZQOLS.fit(
                                y: resampledY,
                                predictors: resampledPredictors,
                                includeConstant: includeConstant,
                                variance: .classical
                            )
                            results.append((replicate, fit.coefficients.map(\.estimate)))
                        }
                        return results
                    }
                }
                var collected = [[Double]](repeating: [], count: reps)
                for try await chunkResults in group {
                    for (replicate, estimates) in chunkResults {
                        collected[replicate] = estimates
                    }
                }
                return collected
            }
        }

        // Singular resamples (NaN rows from the GPU path) are dropped,
        // Stata-style.
        let draws = allDraws.filter { row in !row.contains(where: \.isNaN) }
        let completed = draws.count
        guard completed >= 2 else {
            throw ZQEngineError("bootstrap: too many failed replicates")
        }

        // Bootstrap standard errors: sd of replicate estimates.
        let title = backend == .gpu
            ? "Bootstrap results (pairs, GPU batched)"
            : "Bootstrap results (pairs)"
        var lines = [
            TableFormatter.pad(title, 46, right: false) + "Replications  = " +
            TableFormatter.pad("\(completed)", 10),
            "                                              Number of obs = " +
            TableFormatter.pad("\(n)", 10),
        ]
        if completed < reps {
            lines.append("(\(reps - completed) replicates dropped: singular resample)")
        }
        lines.append("")
        let widths = [12, 12, 12, 8, 8, 22]
        lines.append(
            TableFormatter.pad(sample.responseName, widths[0]) + " | " +
            TableFormatter.pad("Observed", widths[1]) + "  " +
            TableFormatter.pad("Bootstrap", widths[2]) + "  " +
            TableFormatter.pad("z", widths[3]) + "  " +
            TableFormatter.pad("P>|z|", widths[4]) + "  " +
            TableFormatter.pad("[95% conf. interval]", widths[5])
        )
        lines.append(
            TableFormatter.pad("", widths[0]) + " | " +
            TableFormatter.pad("coefficient", widths[1]) + "  " +
            TableFormatter.pad("std. err.", widths[2])
        )
        lines.append(TableFormatter.rule(widths))

        var scalars: [String: Double] = [
            "reps": Double(completed),
            "N": Double(n),
            "gpu": backend == .gpu ? 1 : 0,
        ]
        for j in 0..<k {
            let estimates = draws.map { $0[j] }
            let mean = estimates.reduce(0, +) / Double(completed)
            let variance = estimates.reduce(0) { $0 + ($1 - mean) * ($1 - mean) }
                / Double(completed - 1)
            let se = variance.squareRoot()
            let observed = point.coefficients[j].estimate
            let z = observed / se
            let p = 2 * (1 - ZQDistributions.normalCDF(abs(z)))
            let critical = 1.959963984540054  // Φ⁻¹(0.975)
            lines.append(
                TableFormatter.pad(point.coefficients[j].name, widths[0]) + " | " +
                TableFormatter.pad(TableFormatter.general(observed), widths[1]) + "  " +
                TableFormatter.pad(TableFormatter.general(se), widths[2]) + "  " +
                TableFormatter.pad(TableFormatter.fixed(z, decimals: 2), widths[3]) + "  " +
                TableFormatter.pad(TableFormatter.fixed(p, decimals: 3), widths[4]) + "  " +
                TableFormatter.pad(
                    TableFormatter.general(observed - critical * se) + "   " +
                    TableFormatter.general(observed + critical * se),
                    widths[5]
                )
            )
            scalars["b_\(point.coefficients[j].name)"] = observed
            scalars["se_\(point.coefficients[j].name)"] = se
        }
        return ZQResult(text: lines.joined(separator: "\n"), scalars: scalars)
    }

    // MARK: - Bayesian estimation

    /// `bayes [, mcmcsize(#) burnin(#) seed(#) normalprior(#)]: reg y x…`
    /// — Gibbs-sampled Bayesian linear regression with Stata-style default
    /// priors: N(0, 10000) on coefficients, InvGamma(0.01, 0.01) on σ².
    private func handleBayes(_ command: ZQCommand) throws -> ZQResult {
        guard let body = command.body else {
            throw ZQEngineError("bayes: syntax is 'bayes [, options]: regress …'")
        }
        guard body.verb == "regress" else {
            throw ZQEngineError("bayes currently supports 'regress' bodies only")
        }
        if let seedText = command.option("seed")?.firstArgument {
            guard let value = UInt64(seedText) else {
                throw ZQEngineError("bayes: invalid seed")
            }
            seed = value
        }
        func intOption(_ name: String, default defaultValue: Int) throws -> Int {
            guard let text = command.option(name)?.firstArgument else { return defaultValue }
            guard let value = Int(text), value > 0 else {
                throw ZQEngineError("bayes: invalid \(name)()")
            }
            return value
        }
        let mcmcSize = try intOption("mcmcsize", default: 10_000)
        let burnIn = try intOption("burnin", default: 2_500)
        var priorVariance = 10_000.0
        if let text = command.option("normalprior")?.firstArgument {
            guard let value = Double(text), value > 0 else {
                throw ZQEngineError("bayes: invalid normalprior()")
            }
            priorVariance = value
        }

        let sample = try buildRegressionSample(body, clusterVariable: nil)
        let result = try ZQBayesianRegression.fitGibbs(
            y: sample.y,
            predictors: sample.predictors,
            includeConstant: !body.hasOption("noconstant"),
            mcmcSize: mcmcSize,
            burnIn: burnIn,
            seed: seed,
            coefficientPriorVariance: priorVariance
        )
        lastEstimation = nil   // predict after bayes needs posterior draws

        var lines = ["Bayesian linear regression (Gibbs)"]
        if sample.droppedMissing > 0 {
            lines.append("(\(sample.droppedMissing) observations dropped due to missing values)")
        }
        for (label, value) in [
            ("Number of obs", "\(result.observationCount)"),
            ("MCMC iterations", "\(result.mcmcSize)"),
            ("Burn-in", "\(result.burnIn)"),
            ("Priors", "b ~ N(0, \(TableFormatter.general(priorVariance))), sigma2 ~ IG(.01, .01)"),
        ] {
            lines.append(
                TableFormatter.pad(label, 46, right: false) + "= " +
                TableFormatter.pad(value, 24)
            )
        }
        lines.append("")

        let widths = [12, 12, 12, 24]
        lines.append(
            TableFormatter.pad(sample.responseName, widths[0]) + " | " +
            TableFormatter.pad("Mean", widths[1]) + "  " +
            TableFormatter.pad("Std. dev.", widths[2]) + "  " +
            TableFormatter.pad("[95% cred. interval]", widths[3])
        )
        lines.append(TableFormatter.rule(widths))

        var scalars: [String: Double] = [
            "N": Double(result.observationCount),
            "mcmcsize": Double(result.mcmcSize),
            "burnin": Double(result.burnIn),
        ]
        for coefficient in result.coefficients + [result.sigma] {
            lines.append(
                TableFormatter.pad(coefficient.name, widths[0]) + " | " +
                TableFormatter.pad(TableFormatter.general(coefficient.posteriorMean), widths[1]) + "  " +
                TableFormatter.pad(TableFormatter.general(coefficient.posteriorSD), widths[2]) + "  " +
                TableFormatter.pad(
                    TableFormatter.general(coefficient.credibleLower) + "   " +
                    TableFormatter.general(coefficient.credibleUpper),
                    widths[3]
                )
            )
            scalars["b_\(coefficient.name)"] = coefficient.posteriorMean
            scalars["sd_\(coefficient.name)"] = coefficient.posteriorSD
        }
        return ZQResult(text: lines.joined(separator: "\n"), scalars: scalars)
    }

    // MARK: - Permutation test

    /// `permute, reps(#) [seed(#)]: reg y x…` — permutes the response
    /// within the estimation sample, refits, and reports the empirical
    /// two-sided p-value per coefficient: the fraction of permuted |β*|
    /// at or above the observed |β|. Permutations come from the Philox
    /// argsort stream, so results are exactly reproducible and
    /// replicate-addressable.
    private func handlePermute(_ command: ZQCommand) async throws -> ZQResult {
        guard let body = command.body else {
            throw ZQEngineError("permute: syntax is 'permute, reps(#): command'")
        }
        guard body.verb == "regress" else {
            throw ZQEngineError("permute currently supports 'regress' bodies only")
        }
        guard let repsText = command.option("reps")?.firstArgument,
              let reps = Int(repsText), reps >= 1 else {
            throw ZQEngineError("permute: reps(#) with # ≥ 1 required")
        }
        if let seedText = command.option("seed")?.firstArgument {
            guard let value = UInt64(seedText) else {
                throw ZQEngineError("permute: invalid seed")
            }
            seed = value
        }

        let sample = try buildRegressionSample(body, clusterVariable: nil)
        let includeConstant = !body.hasOption("noconstant")
        let observed = try ZQOLS.fit(
            y: sample.y,
            predictors: sample.predictors,
            includeConstant: includeConstant,
            variance: .classical
        )
        let n = sample.y.count
        let k = observed.coefficients.count
        let generator = Philox4x32(seed: seed)
        let y = sample.y
        let predictors = sample.predictors

        // Count |β*| ≥ |β_obs| per coefficient, in parallel chunks —
        // permutations are counter-addressable so order is irrelevant.
        let observedMagnitudes = observed.coefficients.map { abs($0.estimate) }
        let chunkCount = min(
            max(1, ProcessInfo.processInfo.activeProcessorCount), reps
        )
        let chunkSize = (reps + chunkCount - 1) / chunkCount
        let exceedances: [Int] = try await withThrowingTaskGroup(of: [Int].self) { group in
            for chunk in 0..<chunkCount {
                let lower = chunk * chunkSize
                let upper = min(reps, lower + chunkSize)
                guard lower < upper else { continue }
                group.addTask {
                    var counts = [Int](repeating: 0, count: k)
                    for replicate in lower..<upper {
                        let order = ZQResampling.permutationIndices(
                            replicate: replicate, sampleSize: n, generator: generator
                        )
                        let fit = try ZQOLS.fit(
                            y: order.map { y[$0] },
                            predictors: predictors,
                            includeConstant: includeConstant,
                            variance: .classical
                        )
                        for j in 0..<k
                        where abs(fit.coefficients[j].estimate) >= observedMagnitudes[j] {
                            counts[j] += 1
                        }
                    }
                    return counts
                }
            }
            var total = [Int](repeating: 0, count: k)
            for try await counts in group {
                for j in 0..<k { total[j] += counts[j] }
            }
            return total
        }

        var lines = [
            TableFormatter.pad("Permutation test (response permuted)", 46, right: false) +
            "Replications  = " + TableFormatter.pad("\(reps)", 10),
            "                                              Number of obs = " +
            TableFormatter.pad("\(n)", 10),
            "",
        ]
        let widths = [12, 12, 8, 8, 10, 11]
        lines.append(
            TableFormatter.pad(sample.responseName, widths[0]) + " | " +
            TableFormatter.pad("Observed", widths[1]) + "  " +
            TableFormatter.pad("c", widths[2]) + "  " +
            TableFormatter.pad("reps", widths[3]) + "  " +
            TableFormatter.pad("p=c/reps", widths[4]) + "  " +
            TableFormatter.pad("SE(p)", widths[5])
        )
        lines.append(TableFormatter.rule(widths))

        var scalars: [String: Double] = ["reps": Double(reps), "N": Double(n)]
        for j in 0..<k {
            let name = observed.coefficients[j].name
            // The constant's permutation distribution is degenerate under
            // response permutation; Stata omits it too.
            if name == "_cons" { continue }
            let p = Double(exceedances[j]) / Double(reps)
            let standardError = (p * (1 - p) / Double(reps)).squareRoot()
            lines.append(
                TableFormatter.pad(name, widths[0]) + " | " +
                TableFormatter.pad(
                    TableFormatter.general(observed.coefficients[j].estimate), widths[1]
                ) + "  " +
                TableFormatter.pad("\(exceedances[j])", widths[2]) + "  " +
                TableFormatter.pad("\(reps)", widths[3]) + "  " +
                TableFormatter.pad(TableFormatter.fixed(p, decimals: 4), widths[4]) + "  " +
                TableFormatter.pad(TableFormatter.fixed(standardError, decimals: 4), widths[5])
            )
            scalars["b_\(name)"] = observed.coefficients[j].estimate
            scalars["p_\(name)"] = p
            scalars["c_\(name)"] = Double(exceedances[j])
        }
        return ZQResult(text: lines.joined(separator: "\n"), scalars: scalars)
    }

    // MARK: - Graphics

    private func handleGraph(_ command: ZQCommand) throws -> ZQResult {
        guard !frame.isEmpty else { throw ZQEngineError("no data in memory") }
        let kind: ZQPlotSpec.Kind
        switch command.subverb {
        case "scatter": kind = .scatter
        case "line": kind = .line
        case "histogram": kind = .histogram
        case "kdensity": kind = .kdensity
        case let other:
            throw ZQEngineError("graph \(other ?? ""): not supported yet")
        }

        let names = command.varlist.flatMap(\.referencedNames)
        let mask = try observationMask(command)

        // Single-variable plots: histogram and kernel density.
        if kind == .histogram || kind == .kdensity {
            guard names.count == 1 else {
                throw ZQEngineError("\(command.subverb ?? ""): syntax is '\(command.subverb ?? "") varname'")
            }
            let (values, missing) = try frame.requireNumeric(names[0])
            let data = (0..<frame.rowCount)
                .filter { mask[$0] && !missing[$0] }
                .map { values[$0] }
            guard data.count > 1 else { throw ZQEngineError("no observations") }

            let series: ZQPlotSpec.Series
            if kind == .histogram {
                series = Self.histogramSeries(
                    data: data, label: names[0],
                    binCount: command.option("bins")?.firstArgument.flatMap(Int.init)
                )
            } else {
                series = Self.kernelDensitySeries(data: data, label: names[0])
            }
            lastPlot = ZQPlotSpec(
                kind: kind,
                xLabel: names[0],
                yLabel: kind == .histogram ? "Frequency" : "Density",
                series: [series]
            )
            return ZQResult(
                text: "(plot created: \(data.count) observations)",
                scalars: ["N": Double(data.count)]
            )
        }

        guard names.count == 2 else {
            throw ZQEngineError("graph \(command.subverb ?? ""): syntax is 'graph \(command.subverb ?? "kind") yvar xvar'")
        }
        let (yValues, yMissing) = try frame.requireNumeric(names[0])
        let (xValues, xMissing) = try frame.requireNumeric(names[1])

        func series(label: String, rows: [Int]) -> ZQPlotSpec.Series {
            ZQPlotSpec.Series(
                label: label,
                x: rows.map { xValues[$0] },
                y: rows.map { yValues[$0] }
            )
        }

        let complete = (0..<frame.rowCount).filter {
            mask[$0] && !yMissing[$0] && !xMissing[$0]
        }

        var allSeries: [ZQPlotSpec.Series] = []
        if let byName = command.option("by")?.firstArgument {
            let (byValues, byMissing) = try frame.requireNumeric(byName)
            let levels = ZQFactorExpansion.levels(
                values: complete.map { byValues[$0] },
                missing: complete.map { byMissing[$0] }
            )
            for level in levels {
                let rows = complete.filter { !byMissing[$0] && byValues[$0] == level }
                let label = level == level.rounded()
                    ? "\(byName)=\(Int(level))"
                    : "\(byName)=\(level)"
                allSeries.append(series(label: label, rows: rows))
            }
        } else {
            allSeries.append(series(label: names[0], rows: complete))
        }

        lastPlot = ZQPlotSpec(
            kind: kind,
            xLabel: names[1],
            yLabel: names[0],
            series: allSeries
        )
        let pointCount = allSeries.reduce(0) { $0 + $1.x.count }
        return ZQResult(text: "(plot created: \(pointCount) points)")
    }

    // MARK: - UI support

    /// Compiles filter-bar text into a keep-mask over the current dataset
    /// (the data browser's filter compiles to an `if` expression, §7).
    /// Missing conditions exclude the observation, like `if` qualifiers.
    public func conditionMask(_ expressionText: String) throws -> [Bool] {
        let expression = try parser.parseExpression(expressionText)
        let evaluator = ExpressionEvaluator(frame: frame)
        return try evaluator.evaluateCondition(expression)
    }

    /// Frequency histogram: Sturges bin count by default, half-open bins
    /// with the maximum folded into the last bin. Series points are bin
    /// midpoints vs counts.
    static func histogramSeries(
        data: [Double], label: String, binCount: Int?
    ) -> ZQPlotSpec.Series {
        let n = data.count
        let bins = max(1, binCount ?? Int((Foundation.log2(Double(n))).rounded(.up)) + 1)
        let minimum = data.min()!
        let maximum = data.max()!
        let width = maximum > minimum ? (maximum - minimum) / Double(bins) : 1
        var counts = [Double](repeating: 0, count: bins)
        for value in data {
            let raw = Int((value - minimum) / width)
            counts[min(max(raw, 0), bins - 1)] += 1
        }
        let midpoints = (0..<bins).map { minimum + (Double($0) + 0.5) * width }
        return ZQPlotSpec.Series(label: label, x: midpoints, y: counts)
    }

    /// Epanechnikov kernel density on a 100-point grid with Silverman's
    /// bandwidth h = 0.9·min(sd, IQR/1.349)·n^(−1/5) (Stata's default).
    static func kernelDensitySeries(
        data: [Double], label: String
    ) -> ZQPlotSpec.Series {
        let n = data.count
        let mean = data.reduce(0, +) / Double(n)
        let sd = (data.reduce(0) { $0 + ($1 - mean) * ($1 - mean) }
            / Double(n - 1)).squareRoot()
        let sorted = data.sorted()
        let iqr = ZQSummarize.stataPercentile(sorted: sorted, percent: 75)
            - ZQSummarize.stataPercentile(sorted: sorted, percent: 25)
        var spread = min(sd, iqr / 1.349)
        if spread <= 0 { spread = max(sd, 1e-9) }
        let bandwidth = 0.9 * spread * pow(Double(n), -0.2)

        let lower = sorted.first! - bandwidth
        let upper = sorted.last! + bandwidth
        let gridSize = 100
        let step = (upper - lower) / Double(gridSize - 1)
        var xs = [Double](repeating: 0, count: gridSize)
        var densities = [Double](repeating: 0, count: gridSize)
        for g in 0..<gridSize {
            let point = lower + Double(g) * step
            xs[g] = point
            var total = 0.0
            for value in data {
                let u = (point - value) / bandwidth
                if abs(u) < 1 {
                    total += 0.75 * (1 - u * u)   // Epanechnikov kernel
                }
            }
            densities[g] = total / (Double(n) * bandwidth)
        }
        return ZQPlotSpec.Series(label: label, x: xs, y: densities)
    }

    // MARK: - Qualifiers

    /// Combined if/in keep-mask over current observations.
    private func observationMask(_ command: ZQCommand) throws -> [Bool] {
        var mask = [Bool](repeating: true, count: frame.rowCount)
        if let condition = command.condition {
            let evaluator = ExpressionEvaluator(frame: frame)
            mask = try evaluator.evaluateCondition(condition)
        }
        if let range = command.range {
            let n = frame.rowCount
            func resolve(_ bound: Int) -> Int {
                bound >= 0 ? bound : n + bound + 1
            }
            let lower = max(1, resolve(range.lower))
            let upper = min(n, resolve(range.upper))
            for i in 0..<n where i + 1 < lower || i + 1 > upper {
                mask[i] = false
            }
        }
        return mask
    }
}