spb/metrika Public

1//2//  Session.swift3//  Metrika4//5//  Author:  Simon-Pierre Boucher6//  Contact: contact@spboucher.ai7//  Copyright © 2026 Simon-Pierre Boucher. All rights reserved.8//910import Foundation11import ZQData12import ZQGPU13import ZQGraphics14import ZQParser15import ZQPlanner16import ZQPlugins17import ZQStats1819public struct ZQEngineError: Error, Equatable, Sendable, CustomStringConvertible {20    public var message: String21    public init(_ message: String) { self.message = message }22    public var description: String { message }23}2425/// One interactive session: working dataset, RNG state, panel settings,26/// logging. Every command mutation — console, do-file, or plugin — goes27/// through `execute(_:)` so all frontends share one execution path28/// (CLAUDE.md §7).29public actor ZQSession {30    public private(set) var frame: ZQDataFrame31    public private(set) var lastPlot: ZQPlotSpec?32    public private(set) var seed: UInt64 = 123_456_78933    public private(set) var panelVariable: String?34    public private(set) var timeVariable: String?3536    private let store: ZQDataStore37    private let parser: ZQCommandParser38    private let planner = ZQPlanner(gpuAvailable: ZQGPUBootstrap.isAvailable)39    private let scriptCommands: [String: ZQScriptCommand]40    private let pluginRegistry: ZQPluginRegistry41    private var logFileURL: URL?42    private var scriptDepth = 043    private var lastEstimation: EstimationState?4445    /// - Parameters:46    ///   - discoverUserCommands: scan the commands directory for `.zyq`47    ///     script commands.48    ///   - commandsDirectory: override of the discovery location (tests).49    ///   - plugins: SPM-compiled-in native plugins to register. A plugin50    ///     whose verb shadows a built-in command is rejected.51    public init(52        discoverUserCommands: Bool = true,53        commandsDirectory: URL? = nil,54        plugins: [any ZQCommandPlugin] = []55    ) throws {56        self.store = try ZQDataStore()57        self.frame = try ZQDataFrame()5859        var commands: [String: ZQScriptCommand] = [:]60        if discoverUserCommands {61            let scripts = commandsDirectory.map {62                ZQPluginDiscovery.scriptCommands(in: $0)63            } ?? ZQPluginDiscovery.scriptCommands()64            for script in scripts {65                commands[script.verb] = script66            }67        }68        self.scriptCommands = commands6970        let builtinVerbs = Set(ZQVerbTable.builtin.verbs.keys)71        self.pluginRegistry = try ZQPluginRegistry(72            plugins: plugins, reservedVerbs: builtinVerbs73        )7475        // Plugin verbs join the grammar so their commands parse with full76        // varlist/if/in/options structure (no abbreviations for plugins).77        var verbs = ZQVerbTable.builtin.verbs78        for verb in pluginRegistry.verbs {79            verbs[verb] = verb.count80        }81        self.parser = ZQCommandParser(verbTable: ZQVerbTable(82            verbs: verbs,83            fileVerbs: ZQVerbTable.builtin.fileVerbs,84            assignmentVerbs: ZQVerbTable.builtin.assignmentVerbs,85            prefixVerbs: ZQVerbTable.builtin.prefixVerbs,86            compoundVerbs: ZQVerbTable.builtin.compoundVerbs87        ))88    }8990    // MARK: - Entry point9192    public func execute(_ line: String) async throws -> ZQResult {93        let result: ZQResult94        do {95            result = try await run(line)96        } catch {97            try? appendToLog(command: line, output: "error: \(error)")98            throw error99        }100        try? appendToLog(command: line, output: result.text)101        return result102    }103104    /// Executes a do-file: runs each line, concatenating output. Stops at105    /// the first error, reporting the failing line number.106    public func executeScript(_ text: String) async throws -> ZQResult {107        var outputs: [String] = []108        var scalars: [String: Double] = [:]109        for (index, rawLine) in text.split(110            separator: "\n", omittingEmptySubsequences: false111        ).enumerated() {112            let line = String(rawLine)113            do {114                let result = try await execute(line)115                if !result.text.isEmpty {116                    outputs.append(". \(line)\n\(result.text)")117                }118                scalars.merge(result.scalars) { _, new in new }119            } catch {120                throw ZQEngineError("line \(index + 1): \(error)")121            }122        }123        return ZQResult(text: outputs.joined(separator: "\n\n"), scalars: scalars)124    }125126    private func run(_ line: String) async throws -> ZQResult {127        // Script commands are routed BEFORE the parser: their verbs are128        // not in the grammar and their arguments are raw macro text.129        let trimmed = line.trimmingCharacters(in: .whitespaces)130        let verbWord = String(trimmed.prefix { !$0.isWhitespace })131        if let script = scriptCommands[verbWord] {132            guard scriptDepth < 8 else {133                throw ZQEngineError("user command recursion too deep")134            }135            scriptDepth += 1136            defer { scriptDepth -= 1 }137            let arguments = String(trimmed.dropFirst(verbWord.count))138                .trimmingCharacters(in: .whitespaces)139            let body = script.expandedLines(arguments: arguments)140                .joined(separator: "\n")141            return try await executeScript(body)142        }143144        guard let command = try parser.parse(line) else {145            return ZQResult(text: "")146        }147148        if let registered = pluginRegistry[command.verb] {149            return try await executePlugin(registered, command: command)150        }151152        let plan = planner.plan(command, rowCount: frame.rowCount)153        return try await dispatch(plan.command, backend: plan.backend)154    }155156    /// Runs a native plugin: validates the parsed command against its157    /// declared syntax, hands it a value copy of the dataset, and installs158    /// a returned replacement frame only when `mutates: true`.159    private func executePlugin(160        _ registered: ZQPluginRegistry.Registered, command: ZQCommand161    ) async throws -> ZQResult {162        if let violation = pluginRegistry.validate(command, against: registered.syntax) {163            throw ZQEngineError("\(registered.verb): \(violation)")164        }165        let context = ZQContext(command: command, frame: frame)166        var result = try await registered.plugin.execute(context)167        if let replacement = result.replacementFrame {168            guard registered.syntax.mutates else {169                throw ZQEngineError(170                    "\(registered.verb): plugin returned a dataset but does not declare 'mutates'"171                )172            }173            frame = replacement174            result.replacementFrame = nil175        }176        return result177    }178179    private func dispatch(180        _ command: ZQCommand, backend: ZQBackend = .cpu181    ) async throws -> ZQResult {182        switch command.verb {183        case "use": return try await handleUse(command)184        case "sysuse": return try await handleSysuse(command)185        case "save": return try await handleSave(command)186        case "clear": return handleClear()187        case "describe": return handleDescribe()188        case "summarize": return try handleSummarize(command)189        case "generate": return try handleGenerate(command, replace: false)190        case "replace": return try handleGenerate(command, replace: true)191        case "drop": return try handleDropKeep(command, keep: false)192        case "keep": return try handleDropKeep(command, keep: true)193        case "count": return try handleCount(command)194        case "list": return try handleList(command)195        case "regress": return try handleRegress(command)196        case "xtreg": return try handleXTReg(command)197        case "ivregress": return try handleIVRegress(command)198        case "predict": return try handlePredict(command)199        case "margins": return try handleMargins(command)200        case "elasticnet": return try handleElasticNet(command, defaultAlpha: nil)201        case "boost": return try handleBoost(command)202        case "lasso": return try handleElasticNet(command, defaultAlpha: 1)203        case "logit": return try handleGLM(command, family: .logit)204        case "probit": return try handleGLM(command, family: .probit)205        case "poisson": return try handleGLM(command, family: .poisson)206        case "tabulate": return try handleTabulate(command)207        case "correlate": return try handleCorrelate(command)208        case "set": return try handleSet(command)209        case "xtset": return try handleXTSet(command)210        case "display": return try handleDisplay(command)211        case "bootstrap": return try await handleBootstrap(command, backend: backend)212        case "permute": return try await handlePermute(command)213        case "bayes": return try handleBayes(command)214        case "graph": return try handleGraph(command)215        case "histogram", "scatter", "kdensity":216            var promoted = command217            promoted.subverb = command.verb218            promoted.verb = "graph"219            return try handleGraph(promoted)220        case "log": return try handleLog(command)221        case "help": return handleHelp(command)222        default:223            throw ZQEngineError("command '\(command.verb)' is not implemented yet")224        }225    }226227    // MARK: - Data commands228229    private func handleUse(_ command: ZQCommand) async throws -> ZQResult {230        guard let path = command.argument else {231            throw ZQEngineError("use: file name required")232        }233        let url = URL(fileURLWithPath: (path as NSString).expandingTildeInPath)234        frame = try await store.load(contentsOf: url)235        return ZQResult(236            text: "(\(frame.rowCount) observations, \(frame.columns.count) variables loaded from \(url.lastPathComponent))",237            scalars: ["N": Double(frame.rowCount), "k": Double(frame.columns.count)]238        )239    }240241    /// `sysuse name` — load a bundled sample dataset; bare `sysuse` lists242    /// what ships with the app.243    private func handleSysuse(_ command: ZQCommand) async throws -> ZQResult {244        let samples = ["sales", "mtcars"]245        guard let name = command.argument else {246            return ZQResult(247                text: "Sample datasets: \(samples.joined(separator: ", ")) — load one with 'sysuse <name>'"248            )249        }250        let stem = name.hasSuffix(".csv") ? String(name.dropLast(4)) : name251        guard samples.contains(stem),252              let url = Bundle.module.url(253                forResource: stem, withExtension: "csv", subdirectory: "Samples"254              ) else {255            throw ZQEngineError(256                "sysuse: unknown sample '\(name)' — available: \(samples.joined(separator: ", "))"257            )258        }259        frame = try await store.load(contentsOf: url)260        lastEstimation = nil261        return ZQResult(262            text: "(\(frame.rowCount) observations, \(frame.columns.count) variables — sample dataset '\(stem)')",263            scalars: ["N": Double(frame.rowCount), "k": Double(frame.columns.count)]264        )265    }266267    private func handleSave(_ command: ZQCommand) async throws -> ZQResult {268        guard let path = command.argument else {269            throw ZQEngineError("save: file name required")270        }271        guard !frame.isEmpty else { throw ZQEngineError("no data in memory") }272        let url = URL(fileURLWithPath: (path as NSString).expandingTildeInPath)273        if FileManager.default.fileExists(atPath: url.path),274           !command.hasOption("replace") {275            throw ZQEngineError("file \(url.lastPathComponent) already exists — add ', replace' to overwrite")276        }277        try await store.save(frame, to: url)278        return ZQResult(text: "file \(url.lastPathComponent) saved")279    }280281    private func handleClear() -> ZQResult {282        frame = try! ZQDataFrame()283        panelVariable = nil284        timeVariable = nil285        return ZQResult(text: "")286    }287288    private func handleDescribe() -> ZQResult {289        guard !frame.isEmpty else {290            return ZQResult(text: "Contains data: 0 observations")291        }292        var lines = [293            "Contains data: \(frame.rowCount) observations, \(frame.columns.count) variables",294            "",295            TableFormatter.pad("Variable", 12) + "   " +296            TableFormatter.pad("Type", 8, right: false) + "  Missing",297            String(repeating: "-", count: 34),298        ]299        for column in frame.columns {300            let type = column.data.isNumeric ? "float64" : "string"301            lines.append(302                TableFormatter.pad(column.name, 12) + "   " +303                TableFormatter.pad(type, 8, right: false) + "  " +304                TableFormatter.pad("\(column.data.missingCount)", 7)305            )306        }307        return ZQResult(308            text: lines.joined(separator: "\n"),309            scalars: ["N": Double(frame.rowCount), "k": Double(frame.columns.count)]310        )311    }312313    private func handleSummarize(_ command: ZQCommand) throws -> ZQResult {314        guard !frame.isEmpty else { throw ZQEngineError("no data in memory") }315        let mask = try observationMask(command)316        let detail = command.hasOption("detail")317318        var names: [String] = command.varlist.flatMap(\.referencedNames)319        if names.isEmpty {320            names = frame.columns.filter(\.data.isNumeric).map(\.name)321        }322323        var lines = [324            TableFormatter.pad("Variable", 12) + " | " +325            TableFormatter.pad("Obs", 9) + "  " +326            TableFormatter.pad("Mean", 11) + "  " +327            TableFormatter.pad("Std. dev.", 11) + "  " +328            TableFormatter.pad("Min", 10) + "  " +329            TableFormatter.pad("Max", 10),330            String(repeating: "-", count: 13) + "+" + String(repeating: "-", count: 62),331        ]332        var scalars: [String: Double] = [:]333334        for name in names {335            let (values, missing) = try frame.requireNumeric(name)336            let restrictedMissing = (0..<values.count).map { missing[$0] || !mask[$0] }337            let summary = ZQSummarize.summary(338                name: name, values: values, missing: restrictedMissing, detail: detail339            )340            if let summary {341                lines.append(342                    TableFormatter.pad(name, 12) + " | " +343                    TableFormatter.pad("\(summary.observationCount)", 9) + "  " +344                    TableFormatter.pad(TableFormatter.general(summary.mean), 11) + "  " +345                    TableFormatter.pad(TableFormatter.general(summary.standardDeviation), 11) + "  " +346                    TableFormatter.pad(TableFormatter.general(summary.minimum), 10) + "  " +347                    TableFormatter.pad(TableFormatter.general(summary.maximum), 10)348                )349                if let detailBlock = summary.detail {350                    for p in [1, 5, 10, 25, 50, 75, 90, 95, 99] {351                        lines.append(352                            TableFormatter.pad("p\(p)", 16) + "  " +353                            TableFormatter.pad(354                                TableFormatter.general(detailBlock.percentiles[p] ?? .nan), 11355                            )356                        )357                    }358                    lines.append(359                        TableFormatter.pad("skewness", 16) + "  " +360                        TableFormatter.pad(TableFormatter.general(detailBlock.skewness), 11)361                    )362                    lines.append(363                        TableFormatter.pad("kurtosis", 16) + "  " +364                        TableFormatter.pad(TableFormatter.general(detailBlock.kurtosis), 11)365                    )366                }367                scalars["N"] = Double(summary.observationCount)368                scalars["mean"] = summary.mean369                scalars["sd"] = summary.standardDeviation370                scalars["min"] = summary.minimum371                scalars["max"] = summary.maximum372                scalars["Var"] = summary.variance373            } else {374                lines.append(375                    TableFormatter.pad(name, 12) + " | " +376                    TableFormatter.pad("0", 9)377                )378            }379        }380        return ZQResult(text: lines.joined(separator: "\n"), scalars: scalars)381    }382383    private func handleGenerate(_ command: ZQCommand, replace: Bool) throws -> ZQResult {384        guard let assignment = command.assignment else {385            throw ZQEngineError("\(replace ? "replace" : "generate"): syntax is 'newvar = expression'")386        }387        guard !frame.isEmpty else { throw ZQEngineError("no data in memory") }388389        let evaluator = ExpressionEvaluator(frame: frame)390        var (values, missing) = try evaluator.evaluateNumericColumn(assignment.expression)391        let mask = try observationMask(command)392393        if replace {394            let (oldValues, oldMissing) = try frame.requireNumeric(assignment.target)395            var changed = 0396            for i in 0..<frame.rowCount {397                if mask[i] {398                    let differs = oldMissing[i] != missing[i]399                        || (!missing[i] && oldValues[i] != values[i])400                    if differs { changed += 1 }401                } else {402                    values[i] = oldValues[i]403                    missing[i] = oldMissing[i]404                }405            }406            try frame.replaceColumn(ZQColumn(407                name: assignment.target,408                data: .float64(values: values, missing: missing)409            ))410            return ZQResult(text: "(\(changed) real changes made)")411        } else {412            var missingCount = 0413            for i in 0..<frame.rowCount where !mask[i] {414                values[i] = .nan415                missing[i] = true416            }417            for flag in missing where flag { missingCount += 1 }418            try frame.addColumn(ZQColumn(419                name: assignment.target,420                data: .float64(values: values, missing: missing)421            ))422            let note = missingCount > 0 ? "(\(missingCount) missing values generated)" : ""423            return ZQResult(text: note)424        }425    }426427    private func handleDropKeep(_ command: ZQCommand, keep: Bool) throws -> ZQResult {428        guard !frame.isEmpty else { throw ZQEngineError("no data in memory") }429430        // Observation form: `drop if expr` / `keep if expr` / `keep in 1/10`.431        if command.varlist.isEmpty, command.condition != nil || command.range != nil {432            let mask = try observationMask(command)433            let keepMask = keep ? mask : mask.map { !$0 }434            let before = frame.rowCount435            frame = try frame.filtered(by: keepMask)436            return ZQResult(text: "(\(before - frame.rowCount) observations deleted)")437        }438439        // Variable form: `drop x y` / `keep x y`.440        let names = command.varlist.flatMap(\.referencedNames)441        guard !names.isEmpty else {442            throw ZQEngineError("\(keep ? "keep" : "drop"): specify variables or an 'if' condition")443        }444        if keep {445            let dropped = frame.columnNames.filter { !names.contains($0) }446            try frame.dropColumns(dropped)447        } else {448            try frame.dropColumns(names)449        }450        return ZQResult(text: "")451    }452453    private func handleCount(_ command: ZQCommand) throws -> ZQResult {454        guard !frame.isEmpty else { throw ZQEngineError("no data in memory") }455        let mask = try observationMask(command)456        let count = mask.count { $0 }457        return ZQResult(text: "  \(count)", scalars: ["N": Double(count)])458    }459460    private func handleList(_ command: ZQCommand) throws -> ZQResult {461        guard !frame.isEmpty else { throw ZQEngineError("no data in memory") }462        let mask = try observationMask(command)463        var names = command.varlist.flatMap(\.referencedNames)464        if names.isEmpty { names = frame.columnNames }465        let columns = try names.map { try frame.requireColumn($0) }466467        let rows = (0..<frame.rowCount).filter { mask[$0] }.prefix(50)468        var lines = [names.map { TableFormatter.pad($0, 12) }.joined(separator: "  ")]469        for row in rows {470            let cells = columns.map { column -> String in471                switch column.data {472                case .float64(let values, let missing):473                    return TableFormatter.pad(474                        missing[row] ? "." : TableFormatter.general(values[row]), 12475                    )476                case .string(let values):477                    return TableFormatter.pad(values[row] ?? "", 12)478                }479            }480            lines.append(cells.joined(separator: "  "))481        }482        return ZQResult(text: lines.joined(separator: "\n"))483    }484485    // MARK: - Settings486487    private func handleSet(_ command: ZQCommand) throws -> ZQResult {488        switch command.subverb {489        case "seed":490            guard let text = command.argument, let value = UInt64(text) else {491                throw ZQEngineError("set seed: integer seed required")492            }493            seed = value494            return ZQResult(text: "")495        default:496            throw ZQEngineError("set: unknown setting '\(command.subverb ?? "")'")497        }498    }499500    private func handleXTSet(_ command: ZQCommand) throws -> ZQResult {501        let names = command.varlist.flatMap(\.referencedNames)502        guard (1...2).contains(names.count) else {503            throw ZQEngineError("xtset: syntax is 'xtset panelvar [timevar]'")504        }505        for name in names {506            _ = try frame.requireColumn(name)507        }508        panelVariable = names[0]509        timeVariable = names.count > 1 ? names[1] : nil510        var text = "Panel variable: \(names[0])"511        if let timeVariable { text += "\nTime variable: \(timeVariable)" }512        return ZQResult(text: text)513    }514515    private func handleDisplay(_ command: ZQCommand) throws -> ZQResult {516        guard let assignment = command.assignment else {517            throw ZQEngineError("display: expression required")518        }519        let evaluator = ExpressionEvaluator(frame: frame)520        switch try evaluator.evaluate(assignment.expression) {521        case .numericScalar(let value):522            return ZQResult(text: TableFormatter.general(value, significant: 10))523        case .stringScalar(let value):524            return ZQResult(text: value)525        case .missingScalar:526            return ZQResult(text: ".")527        case .numeric(let values, let missing):528            guard let first = values.first else { return ZQResult(text: ".") }529            return ZQResult(530                text: missing[0] ? "." : TableFormatter.general(first, significant: 10)531            )532        case .strings(let values):533            return ZQResult(text: values.first.flatMap { $0 } ?? "")534        }535    }536537    private func handleLog(_ command: ZQCommand) throws -> ZQResult {538        if command.argument == "close" || command.subverb == "close" {539            logFileURL = nil540            return ZQResult(text: "(log closed)")541        }542        guard var path = command.argument else {543            throw ZQEngineError("log: syntax is 'log using filename' or 'log close'")544        }545        if path == "close" {546            logFileURL = nil547            return ZQResult(text: "(log closed)")548        }549        if path.hasPrefix("using ") { path = String(path.dropFirst(6)) }550        let url = URL(fileURLWithPath: (path as NSString).expandingTildeInPath)551        if !FileManager.default.fileExists(atPath: url.path) {552            FileManager.default.createFile(atPath: url.path, contents: nil)553        }554        logFileURL = url555        return ZQResult(text: "(log started: \(url.path))")556    }557558    /// `help [command]` — renders the shared command reference.559    private func handleHelp(_ command: ZQCommand) -> ZQResult {560        let requested = command.varlist.flatMap(\.referencedNames).first561            ?? command.argument562563        if let requested {564            guard let doc = ZQCommandReference.doc(for: requested) else {565                return ZQResult(text: "help: no entry for '\(requested)' — type 'help' for the full list")566            }567            var lines = [568                "\(doc.verb)\(doc.abbreviation.map { " (abbreviation: \($0))" } ?? "")",569                String(repeating: "-", count: 60),570                doc.summary,571                "",572                "Syntax:  \(doc.syntax)",573            ]574            if !doc.options.isEmpty {575                lines.append("")576                lines.append("Options:")577                for option in doc.options {578                    lines.append("  " + TableFormatter.pad(option.name, 16, right: false) + option.meaning)579                }580            }581            if !doc.examples.isEmpty {582                lines.append("")583                lines.append("Examples:")584                for example in doc.examples { lines.append("  . \(example)") }585            }586            if let notes = doc.notes {587                lines.append("")588                lines.append(notes)589            }590            return ZQResult(text: lines.joined(separator: "\n"))591        }592593        var lines = ["Metrika command reference — 'help <command>' for details", ""]594        for category in ZQCommandReference.categories {595            let docs = ZQCommandReference.all.filter { $0.category == category }596            guard !docs.isEmpty else { continue }597            lines.append(category)598            for doc in docs {599                lines.append(600                    "  " + TableFormatter.pad(doc.verb, 14, right: false) + doc.summary.prefix(60)601                )602            }603            lines.append("")604        }605        return ZQResult(text: lines.joined(separator: "\n"))606    }607608    private func appendToLog(command: String, output: String) throws {609        guard let logFileURL else { return }610        let entry = ". \(command)\n\(output)\n\n"611        let handle = try FileHandle(forWritingTo: logFileURL)612        defer { try? handle.close() }613        try handle.seekToEnd()614        try handle.write(contentsOf: Data(entry.utf8))615    }616617    // MARK: - Estimation618619    /// How a fitted regressor is recomputed from the working dataset —620    /// the estimation-state piece that lets `predict` (and later621    /// `margins`) evaluate xb on arbitrary observations.622    enum PredictorDefinition: Equatable, Sendable {623        case column(String)624        case indicator(variable: String, level: Double)625        case product([String])626        case constant627    }628629    /// Stata's e() analog: what the last estimation command fitted.630    struct EstimationState: Sendable {631        enum Kind: Equatable, Sendable {632            case ols, iv633            case glm(ZQGLMFamily)634            case boost(ZQBoostModel)635        }636637        var kind: Kind638        var responseName: String639        /// Coefficients aligned with their recomputation recipes.640        var coefficients: [(name: String, value: Double, definition: PredictorDefinition)]641        /// Covariance matrix, column-major k×k aligned with `coefficients`.642        var vce: [Double]643        /// Degrees of freedom for t inference; nil means normal (ML).644        var inferenceDF: Double?645        /// Estimation-sample design (column-major n×k incl. constant) —646        /// kept for GLMs, whose average marginal effects need it.647        var sampleDesign: [Double]?648        var sampleSize: Int649    }650651    /// Assembled estimation sample after if/in restriction and listwise652    /// deletion (§6: explicit report line for dropped observations).653    struct RegressionSample {654        var responseName: String655        var y: [Double]656        var predictors: [(name: String, values: [Double])]657        /// Aligned with `predictors`.658        var definitions: [PredictorDefinition]659        var clusterLabels: [Int]?660        var droppedMissing: Int661    }662663    private func buildRegressionSample(664        _ command: ZQCommand, clusterVariable: String?665    ) throws -> RegressionSample {666        guard !frame.isEmpty else { throw ZQEngineError("no data in memory") }667        guard let first = command.varlist.first else {668            throw ZQEngineError("\(command.verb): dependent variable required")669        }670        guard case .simple(let responseName) = first else {671            throw ZQEngineError("\(command.verb): dependent variable cannot use factor notation")672        }673674        let mask = try observationMask(command)675        let (yAll, yMissing) = try frame.requireNumeric(responseName)676677        // Collect raw regressor sources (numeric columns) for missingness.678        struct Source {679            var spec: ZQVarSpec680            var columns: [(name: String, values: [Double], missing: [Bool])]681        }682        var sources: [Source] = []683        for spec in command.varlist.dropFirst() {684            var columns: [(String, [Double], [Bool])] = []685            for name in spec.referencedNames {686                let (values, missing) = try frame.requireNumeric(name)687                columns.append((name, values, missing))688            }689            sources.append(Source(spec: spec, columns: columns))690        }691692        var clusterColumn: ZQColumn?693        if let clusterVariable {694            clusterColumn = try frame.requireColumn(clusterVariable)695        }696697        // Listwise deletion over response, regressors, and cluster var.698        var keptRows: [Int] = []699        var droppedMissing = 0700        for i in 0..<frame.rowCount where mask[i] {701            var complete = !yMissing[i]702            if complete {703                outer: for source in sources {704                    for column in source.columns where column.missing[i] {705                        complete = false706                        break outer707                    }708                }709            }710            if complete, let clusterColumn {711                switch clusterColumn.data {712                case .float64(_, let missing): complete = !missing[i]713                case .string(let values): complete = values[i] != nil714                }715            }716            if complete {717                keptRows.append(i)718            } else {719                droppedMissing += 1720            }721        }722        guard !keptRows.isEmpty else {723            throw ZQEngineError("no observations")724        }725726        let y = keptRows.map { yAll[$0] }727728        // Expand specs on the estimation sample (plan-time expansion),729        // recording each column's recomputation recipe for predict.730        var predictors: [(name: String, values: [Double])] = []731        var definitions: [PredictorDefinition] = []732        for source in sources {733            switch source.spec {734            case .simple(let name):735                let (values, _) = try frame.requireNumeric(name)736                predictors.append((name, keptRows.map { values[$0] }))737                definitions.append(.column(name))738            case .factor(let op, let name):739                let (values, missing) = try frame.requireNumeric(name)740                switch op {741                case "i", "ib":742                    let sampleValues = keptRows.map { values[$0] }743                    let sampleMissing = keptRows.map { missing[$0] }744                    let expanded = ZQFactorExpansion.expand(745                        name: name, values: sampleValues, missing: sampleMissing746                    )747                    guard !expanded.isEmpty else {748                        throw ZQEngineError("factor variable '\(name)' has a single level")749                    }750                    for term in expanded {751                        predictors.append((term.name, term.values))752                        definitions.append(.indicator(variable: name, level: term.level))753                    }754                case "c":755                    predictors.append((name, keptRows.map { values[$0] }))756                    definitions.append(.column(name))757                default:758                    throw ZQEngineError("unsupported factor operator '\(op).'")759                }760            case .interaction(let parts, _):761                // Continuous-by-continuous interactions only.762                var label: [String] = []763                var factors: [String] = []764                var product = [Double](repeating: 1, count: keptRows.count)765                for part in parts {766                    guard case .factor(let op, let name) = part, op == "c" else {767                        guard case .simple(let name) = part else {768                            throw ZQEngineError("factor interactions beyond c.#c. are not supported yet")769                        }770                        let (values, _) = try frame.requireNumeric(name)771                        for (j, row) in keptRows.enumerated() { product[j] *= values[row] }772                        label.append(name)773                        factors.append(name)774                        continue775                    }776                    let (values, _) = try frame.requireNumeric(name)777                    for (j, row) in keptRows.enumerated() { product[j] *= values[row] }778                    label.append("c.\(name)")779                    factors.append(name)780                }781                predictors.append((label.joined(separator: "#"), product))782                definitions.append(.product(factors))783            }784        }785786        // Cluster labels → dense integer codes.787        var clusterLabels: [Int]?788        if let clusterColumn {789            var codes: [Int] = []790            codes.reserveCapacity(keptRows.count)791            switch clusterColumn.data {792            case .float64(let values, _):793                var mapping: [Double: Int] = [:]794                for row in keptRows {795                    let value = values[row]796                    if let code = mapping[value] {797                        codes.append(code)798                    } else {799                        mapping[value] = mapping.count800                        codes.append(mapping.count - 1)801                    }802                }803            case .string(let values):804                var mapping: [String: Int] = [:]805                for row in keptRows {806                    let value = values[row] ?? ""807                    if let code = mapping[value] {808                        codes.append(code)809                    } else {810                        mapping[value] = mapping.count811                        codes.append(mapping.count - 1)812                    }813                }814            }815            clusterLabels = codes816        }817818        return RegressionSample(819            responseName: responseName,820            y: y,821            predictors: predictors,822            definitions: definitions,823            clusterLabels: clusterLabels,824            droppedMissing: droppedMissing825        )826    }827828    /// Records the estimation state (Stata's e()) for post-estimation829    /// commands. The constant, when fitted, is the last coefficient.830    private func recordEstimation(831        kind: EstimationState.Kind,832        responseName: String,833        coefficients: [ZQCoefficient],834        definitions: [PredictorDefinition],835        includeConstant: Bool,836        vce: [Double] = [],837        inferenceDF: Double? = nil,838        sampleDesign: [Double]? = nil,839        sampleSize: Int = 0840    ) {841        var recipes = definitions842        if includeConstant { recipes.append(.constant) }843        guard recipes.count == coefficients.count else {844            lastEstimation = nil845            return846        }847        lastEstimation = EstimationState(848            kind: kind,849            responseName: responseName,850            coefficients: zip(coefficients, recipes).map {851                (name: $0.name, value: $0.estimate, definition: $1)852            },853            vce: vce,854            inferenceDF: inferenceDF,855            sampleDesign: sampleDesign,856            sampleSize: sampleSize857        )858    }859860    /// Column-major estimation-sample design from sample predictors.861    private func designMatrix(862        from sample: RegressionSample, includeConstant: Bool863    ) -> [Double] {864        let n = sample.y.count865        var design = [Double]()866        design.reserveCapacity(n * (sample.predictors.count + 1))867        for column in sample.predictors { design.append(contentsOf: column.values) }868        if includeConstant { design.append(contentsOf: [Double](repeating: 1, count: n)) }869        return design870    }871872    private func varianceEstimator(873        _ command: ZQCommand, clusterLabels: [Int]?874    ) throws -> ZQVarianceEstimator {875        if let clusterLabels { return .cluster(clusterLabels) }876        if command.hasOption("hc2") { return .hc2 }877        if command.hasOption("hc3") { return .hc3 }878        if command.hasOption("robust") || command.hasOption("r") { return .hc1 }879        if let vce = command.option("vce")?.firstArgument {880            switch vce {881            case "robust": return .hc1882            case "hc2": return .hc2883            case "hc3": return .hc3884            default:885                throw ZQEngineError("vce(\(vce)) is not supported")886            }887        }888        return .classical889    }890891    private func handleRegress(_ command: ZQCommand) throws -> ZQResult {892        let clusterVariable = command.option("cluster")?.firstArgument893            ?? clusterFromVCE(command)894        let sample = try buildRegressionSample(command, clusterVariable: clusterVariable)895        let variance = try varianceEstimator(command, clusterLabels: sample.clusterLabels)896        let level = try confidenceLevel(command)897898        let result = try ZQOLS.fit(899            y: sample.y,900            predictors: sample.predictors,901            includeConstant: !command.hasOption("noconstant"),902            variance: variance,903            confidenceLevel: level904        )905        recordEstimation(906            kind: .ols,907            responseName: sample.responseName,908            coefficients: result.coefficients,909            definitions: sample.definitions,910            includeConstant: !command.hasOption("noconstant"),911            vce: result.vce,912            inferenceDF: result.inferenceDF,913            sampleSize: result.observationCount914        )915916        var header: [String] = []917        switch variance {918        case .classical: header.append("Linear regression")919        case .cluster: header.append("Linear regression (cluster-robust)")920        default: header.append("Linear regression (robust)")921        }922        if sample.droppedMissing > 0 {923            header.append("(\(sample.droppedMissing) observations dropped due to missing values)")924        }925926        var stats = [927            ("Number of obs", TableFormatter.general(Double(result.observationCount))),928            ("R-squared", TableFormatter.fixed(result.rSquared, decimals: 4)),929            ("Adj R-squared", TableFormatter.fixed(result.adjustedRSquared, decimals: 4)),930            ("Root MSE", TableFormatter.general(result.rootMSE, significant: 6)),931        ]932        if let f = result.fStatistic, let df = result.fDF, let p = result.fPValue {933            stats.insert(934                ("F(\(Int(df.0)), \(Int(df.1)))", TableFormatter.fixed(f, decimals: 2)),935                at: 1936            )937            stats.insert(("Prob > F", TableFormatter.fixed(p, decimals: 4)), at: 2)938        }939        if let g = result.clusterCount {940            stats.append(("Clusters", "\(g)"))941        }942        for (label, value) in stats {943            header.append(944                TableFormatter.pad(label, 46, right: false) + "= " +945                TableFormatter.pad(value, 12)946            )947        }948949        var lines = header950        lines.append("")951        let widths = [12, 12, 11, 8, 8, 22]952        lines.append(953            TableFormatter.pad(sample.responseName, widths[0]) + " | " +954            TableFormatter.pad("Coefficient", widths[1]) + "  " +955            TableFormatter.pad("Std. err.", widths[2]) + "  " +956            TableFormatter.pad("t", widths[3]) + "  " +957            TableFormatter.pad("P>|t|", widths[4]) + "  " +958            TableFormatter.pad("[\(Int(level * 100))% conf. interval]", widths[5])959        )960        lines.append(TableFormatter.rule(widths))961        var scalars: [String: Double] = [962            "N": Double(result.observationCount),963            "r2": result.rSquared,964            "r2_a": result.adjustedRSquared,965            "rmse": result.rootMSE,966            "df_r": result.inferenceDF,967        ]968        if let f = result.fStatistic { scalars["F"] = f }969        if let g = result.clusterCount { scalars["N_clust"] = Double(g) }970971        for coefficient in result.coefficients {972            lines.append(973                TableFormatter.pad(coefficient.name, widths[0]) + " | " +974                TableFormatter.pad(TableFormatter.general(coefficient.estimate), widths[1]) + "  " +975                TableFormatter.pad(TableFormatter.general(coefficient.standardError), widths[2]) + "  " +976                TableFormatter.pad(TableFormatter.fixed(coefficient.tStatistic, decimals: 2), widths[3]) + "  " +977                TableFormatter.pad(TableFormatter.fixed(coefficient.pValue, decimals: 3), widths[4]) + "  " +978                TableFormatter.pad(979                    TableFormatter.general(coefficient.confidenceLower) + "   " +980                    TableFormatter.general(coefficient.confidenceUpper),981                    widths[5]982                )983            )984            scalars["b_\(coefficient.name)"] = coefficient.estimate985            scalars["se_\(coefficient.name)"] = coefficient.standardError986        }987        return ZQResult(text: lines.joined(separator: "\n"), scalars: scalars)988    }989990    private func handleXTReg(_ command: ZQCommand) throws -> ZQResult {991        guard command.hasOption("fe") else {992            throw ZQEngineError("xtreg: only the fixed-effects estimator is implemented — add ', fe'")993        }994        guard let panelVariable else {995            throw ZQEngineError("xtreg: declare the panel first with 'xtset panelvar [timevar]'")996        }997        // v0.2 restriction: the cluster variable must be the panel variable.998        var clustered = false999        if let clusterName = command.option("cluster")?.firstArgument1000            ?? clusterFromVCE(command) {1001            guard clusterName == panelVariable else {1002                throw ZQEngineError(1003                    "xtreg: clustering on a variable other than the panel (\(panelVariable)) is not supported yet"1004                )1005            }1006            clustered = true1007        }10081009        // Panel codes ride the same listwise-deletion path as cluster codes.1010        let sample = try buildRegressionSample(command, clusterVariable: panelVariable)1011        guard let panelCodes = sample.clusterLabels else {1012            throw ZQEngineError("xtreg: panel variable missing from sample")1013        }10141015        let result = try ZQFixedEffects.fitWithin(1016            y: sample.y,1017            predictors: sample.predictors,1018            groups: panelCodes,1019            clustered: clustered,1020            confidenceLevel: try confidenceLevel(command)1021        )1022        // predict after xtreg needs the estimated u_i — not stored yet.1023        lastEstimation = nil10241025        var header = ["Fixed-effects (within) regression"]1026        if sample.droppedMissing > 0 {1027            header.append("(\(sample.droppedMissing) observations dropped due to missing values)")1028        }1029        var stats: [(String, String)] = [1030            ("Number of obs", "\(result.observationCount)"),1031            ("Number of groups", "\(result.groupCount)"),1032            ("R-squared (within)", TableFormatter.fixed(result.rSquaredWithin, decimals: 4)),1033            ("Root MSE", TableFormatter.general(result.rootMSE, significant: 6)),1034        ]1035        if let f = result.fStatistic, let p = result.fPValue {1036            stats.append(("F", TableFormatter.fixed(f, decimals: 2)))1037            stats.append(("Prob > F", TableFormatter.fixed(p, decimals: 4)))1038        }1039        for (label, value) in stats {1040            header.append(1041                TableFormatter.pad(label, 46, right: false) + "= " +1042                TableFormatter.pad(value, 12)1043            )1044        }10451046        var scalars: [String: Double] = [1047            "N": Double(result.observationCount),1048            "N_g": Double(result.groupCount),1049            "r2_w": result.rSquaredWithin,1050            "rmse": result.rootMSE,1051            "df_r": result.inferenceDF,1052        ]1053        let lines = header + [""] + coefficientTable(1054            response: sample.responseName,1055            coefficients: result.coefficients,1056            statisticLabel: "t",1057            level: try confidenceLevel(command),1058            scalars: &scalars1059        )1060        return ZQResult(text: lines.joined(separator: "\n"), scalars: scalars)1061    }10621063    private func handleIVRegress(_ command: ZQCommand) throws -> ZQResult {1064        guard command.subverb == "2sls" else {1065            throw ZQEngineError("ivregress: only the 2sls estimator is implemented")1066        }1067        guard let ivSpec = command.ivSpec else {1068            throw ZQEngineError("ivregress: syntax is 'ivregress 2sls depvar [exog] (endog = instruments)'")1069        }1070        guard !frame.isEmpty else { throw ZQEngineError("no data in memory") }1071        guard let first = command.varlist.first, case .simple(let response) = first else {1072            throw ZQEngineError("ivregress: dependent variable required")1073        }1074        let exogenousNames = command.varlist.dropFirst().flatMap(\.referencedNames)1075        let clusterVariable = command.option("cluster")?.firstArgument1076            ?? clusterFromVCE(command)10771078        // Listwise deletion across every variable involved.1079        let allNames = [response] + exogenousNames + ivSpec.endogenous + ivSpec.instruments1080        let sources = try allNames.map { try frame.requireNumeric($0) }1081        var clusterColumn: ZQColumn?1082        if let clusterVariable {1083            clusterColumn = try frame.requireColumn(clusterVariable)1084        }1085        let mask = try observationMask(command)1086        var keptRows: [Int] = []1087        var dropped = 01088        for i in 0..<frame.rowCount where mask[i] {1089            var complete = sources.allSatisfy { !$0.missing[i] }1090            if complete, let clusterColumn {1091                switch clusterColumn.data {1092                case .float64(_, let missing): complete = !missing[i]1093                case .string(let values): complete = values[i] != nil1094                }1095            }1096            if complete { keptRows.append(i) } else { dropped += 1 }1097        }1098        guard !keptRows.isEmpty else { throw ZQEngineError("no observations") }10991100        func gather(_ names: [String]) throws -> [(name: String, values: [Double])] {1101            try names.map { name in1102                let (values, _) = try frame.requireNumeric(name)1103                return (name, keptRows.map { values[$0] })1104            }1105        }1106        let y = try frame.requireNumeric(response).values1107        let clusterLabels: [Int]? = try clusterVariable.map { name in1108            let column = try frame.requireColumn(name)1109            var mapping: [String: Int] = [:]1110            var codes: [Int] = []1111            for row in keptRows {1112                let key: String1113                switch column.data {1114                case .float64(let values, _): key = "\(values[row])"1115                case .string(let values): key = values[row] ?? ""1116                }1117                if let code = mapping[key] {1118                    codes.append(code)1119                } else {1120                    mapping[key] = mapping.count1121                    codes.append(mapping.count - 1)1122                }1123            }1124            return codes1125        }11261127        let result = try ZQIV.fit2SLS(1128            y: keptRows.map { y[$0] },1129            endogenous: try gather(ivSpec.endogenous),1130            exogenous: try gather(exogenousNames),1131            instruments: try gather(ivSpec.instruments),1132            includeConstant: !command.hasOption("noconstant"),1133            variance: try varianceEstimator(command, clusterLabels: clusterLabels),1134            confidenceLevel: try confidenceLevel(command)1135        )1136        recordEstimation(1137            kind: .iv,1138            responseName: response,1139            coefficients: result.coefficients,1140            definitions: (ivSpec.endogenous + exogenousNames).map { .column($0) },1141            includeConstant: !command.hasOption("noconstant"),1142            vce: result.vce,1143            inferenceDF: result.inferenceDF,1144            sampleSize: result.observationCount1145        )11461147        var header = ["Instrumental variables (2SLS) regression"]1148        if dropped > 0 {1149            header.append("(\(dropped) observations dropped due to missing values)")1150        }1151        header.append("Instrumented: \(ivSpec.endogenous.joined(separator: " "))")1152        header.append("Instruments:  \(ivSpec.instruments.joined(separator: " "))")1153        var stats: [(String, String)] = [1154            ("Number of obs", "\(result.observationCount)"),1155            ("R-squared", TableFormatter.fixed(result.rSquared, decimals: 4)),1156            ("Root MSE", TableFormatter.general(result.rootMSE, significant: 6)),1157        ]1158        if let f = result.fStatistic, let p = result.fPValue {1159            stats.insert(("F", TableFormatter.fixed(f, decimals: 2)), at: 1)1160            stats.insert(("Prob > F", TableFormatter.fixed(p, decimals: 4)), at: 2)1161        }1162        if let g = result.clusterCount { stats.append(("Clusters", "\(g)")) }1163        for (label, value) in stats {1164            header.append(1165                TableFormatter.pad(label, 46, right: false) + "= " +1166                TableFormatter.pad(value, 12)1167            )1168        }11691170        var scalars: [String: Double] = [1171            "N": Double(result.observationCount),1172            "r2": result.rSquared,1173            "rmse": result.rootMSE,1174            "df_r": result.inferenceDF,1175        ]1176        if let f = result.fStatistic { scalars["F"] = f }1177        let lines = header + [""] + coefficientTable(1178            response: response,1179            coefficients: result.coefficients,1180            statisticLabel: "t",1181            level: try confidenceLevel(command),1182            scalars: &scalars1183        )1184        return ZQResult(text: lines.joined(separator: "\n"), scalars: scalars)1185    }11861187    /// Shared coefficient-table rendering (regress/xtreg/ivregress/GLMs).1188    private func coefficientTable(1189        response: String,1190        coefficients: [ZQCoefficient],1191        statisticLabel: String,1192        level: Double,1193        scalars: inout [String: Double]1194    ) -> [String] {1195        let widths = [12, 12, 11, 8, 8, 22]1196        var lines = [1197            TableFormatter.pad(response, widths[0]) + " | " +1198            TableFormatter.pad("Coefficient", widths[1]) + "  " +1199            TableFormatter.pad("Std. err.", widths[2]) + "  " +1200            TableFormatter.pad(statisticLabel, widths[3]) + "  " +1201            TableFormatter.pad("P>|\(statisticLabel)|", widths[4]) + "  " +1202            TableFormatter.pad("[\(Int(level * 100))% conf. interval]", widths[5]),1203            TableFormatter.rule(widths),1204        ]1205        for coefficient in coefficients {1206            lines.append(1207                TableFormatter.pad(coefficient.name, widths[0]) + " | " +1208                TableFormatter.pad(TableFormatter.general(coefficient.estimate), widths[1]) + "  " +1209                TableFormatter.pad(TableFormatter.general(coefficient.standardError), widths[2]) + "  " +1210                TableFormatter.pad(TableFormatter.fixed(coefficient.tStatistic, decimals: 2), widths[3]) + "  " +1211                TableFormatter.pad(TableFormatter.fixed(coefficient.pValue, decimals: 3), widths[4]) + "  " +1212                TableFormatter.pad(1213                    TableFormatter.general(coefficient.confidenceLower) + "   " +1214                    TableFormatter.general(coefficient.confidenceUpper),1215                    widths[5]1216                )1217            )1218            scalars["b_\(coefficient.name)"] = coefficient.estimate1219            scalars["se_\(coefficient.name)"] = coefficient.standardError1220        }1221        return lines1222    }12231224    private func handleGLM(1225        _ command: ZQCommand, family: ZQGLMFamily1226    ) throws -> ZQResult {1227        let clusterVariable = command.option("cluster")?.firstArgument1228            ?? clusterFromVCE(command)1229        let sample = try buildRegressionSample(command, clusterVariable: clusterVariable)1230        let variance = try varianceEstimator(command, clusterLabels: sample.clusterLabels)1231        let level = try confidenceLevel(command)12321233        let result = try ZQGLM.fit(1234            y: sample.y,1235            predictors: sample.predictors,1236            family: family,1237            includeConstant: !command.hasOption("noconstant"),1238            variance: variance,1239            confidenceLevel: level1240        )1241        recordEstimation(1242            kind: .glm(family),1243            responseName: sample.responseName,1244            coefficients: result.coefficients,1245            definitions: sample.definitions,1246            includeConstant: !command.hasOption("noconstant"),1247            vce: result.vce,1248            sampleDesign: designMatrix(1249                from: sample,1250                includeConstant: !command.hasOption("noconstant")1251            ),1252            sampleSize: result.observationCount1253        )12541255        let title: String1256        switch family {1257        case .logit: title = "Logistic regression"1258        case .probit: title = "Probit regression"1259        case .poisson: title = "Poisson regression"1260        }12611262        var header = [title]1263        if sample.droppedMissing > 0 {1264            header.append("(\(sample.droppedMissing) observations dropped due to missing values)")1265        }1266        var stats: [(String, String)] = [1267            ("Number of obs", "\(result.observationCount)"),1268        ]1269        if let chi2 = result.chiSquared, let p = result.chiSquaredPValue {1270            let label: String1271            if case .classical = variance {1272                label = "LR chi2(\(result.chiSquaredDF))"1273            } else {1274                label = "Wald chi2(\(result.chiSquaredDF))"1275            }1276            stats.append((label, TableFormatter.fixed(chi2, decimals: 2)))1277            stats.append(("Prob > chi2", TableFormatter.fixed(p, decimals: 4)))1278        }1279        stats.append(("Log likelihood", TableFormatter.fixed(result.logLikelihood, decimals: 4)))1280        stats.append(("Pseudo R2", TableFormatter.fixed(result.pseudoRSquared, decimals: 4)))1281        if let g = result.clusterCount {1282            stats.append(("Clusters", "\(g)"))1283        }1284        for (label, value) in stats {1285            header.append(1286                TableFormatter.pad(label, 46, right: false) + "= " +1287                TableFormatter.pad(value, 12)1288            )1289        }12901291        var lines = header1292        lines.append("")1293        let widths = [12, 12, 11, 8, 8, 22]1294        lines.append(1295            TableFormatter.pad(sample.responseName, widths[0]) + " | " +1296            TableFormatter.pad("Coefficient", widths[1]) + "  " +1297            TableFormatter.pad("Std. err.", widths[2]) + "  " +1298            TableFormatter.pad("z", widths[3]) + "  " +1299            TableFormatter.pad("P>|z|", widths[4]) + "  " +1300            TableFormatter.pad("[\(Int(level * 100))% conf. interval]", widths[5])1301        )1302        lines.append(TableFormatter.rule(widths))13031304        var scalars: [String: Double] = [1305            "N": Double(result.observationCount),1306            "ll": result.logLikelihood,1307            "ll_0": result.nullLogLikelihood,1308            "r2_p": result.pseudoRSquared,1309        ]1310        if let chi2 = result.chiSquared { scalars["chi2"] = chi2 }1311        if let g = result.clusterCount { scalars["N_clust"] = Double(g) }13121313        for coefficient in result.coefficients {1314            lines.append(1315                TableFormatter.pad(coefficient.name, widths[0]) + " | " +1316                TableFormatter.pad(TableFormatter.general(coefficient.estimate), widths[1]) + "  " +1317                TableFormatter.pad(TableFormatter.general(coefficient.standardError), widths[2]) + "  " +1318                TableFormatter.pad(TableFormatter.fixed(coefficient.tStatistic, decimals: 2), widths[3]) + "  " +1319                TableFormatter.pad(TableFormatter.fixed(coefficient.pValue, decimals: 3), widths[4]) + "  " +1320                TableFormatter.pad(1321                    TableFormatter.general(coefficient.confidenceLower) + "   " +1322                    TableFormatter.general(coefficient.confidenceUpper),1323                    widths[5]1324                )1325            )1326            scalars["b_\(coefficient.name)"] = coefficient.estimate1327            scalars["se_\(coefficient.name)"] = coefficient.standardError1328        }1329        return ZQResult(text: lines.joined(separator: "\n"), scalars: scalars)1330    }13311332    // MARK: - Tabulate & correlate13331334    private func handleTabulate(_ command: ZQCommand) throws -> ZQResult {1335        guard !frame.isEmpty else { throw ZQEngineError("no data in memory") }1336        let names = command.varlist.flatMap(\.referencedNames)1337        guard names.count == 1 || names.count == 2 else {1338            throw ZQEngineError("tabulate: syntax is 'tabulate varname [varname2]'")1339        }1340        let mask = try observationMask(command)13411342        /// Rendered level label per observation; nil = missing (excluded1343        /// unless the `missing` option is given).1344        func labels(_ name: String) throws -> [String?] {1345            let column = try frame.requireColumn(name)1346            switch column.data {1347            case .float64(let values, let missing):1348                return (0..<frame.rowCount).map { i in1349                    missing[i] ? nil : TableFormatter.general(values[i])1350                }1351            case .string(let values):1352                return values1353            }1354        }13551356        let includeMissing = command.hasOption("missing")1357        let first = try labels(names[0])13581359        if names.count == 1 {1360            var counts: [String: Int] = [:]1361            var total = 01362            for i in 0..<frame.rowCount where mask[i] {1363                guard let label = includeMissing ? (first[i] ?? ".") : first[i] else {1364                    continue1365                }1366                counts[label, default: 0] += 11367                total += 11368            }1369            guard total > 0 else { throw ZQEngineError("no observations") }13701371            var lines = [1372                TableFormatter.pad(names[0], 14) + " | " +1373                TableFormatter.pad("Freq.", 9) + "  " +1374                TableFormatter.pad("Percent", 9) + "  " +1375                TableFormatter.pad("Cum.", 9),1376                String(repeating: "-", count: 15) + "+" + String(repeating: "-", count: 33),1377            ]1378            var cumulative = 0.01379            for label in counts.keys.sorted(by: numericAwareLess) {1380                let count = counts[label]!1381                let percent = 100 * Double(count) / Double(total)1382                cumulative += percent1383                lines.append(1384                    TableFormatter.pad(label, 14) + " | " +1385                    TableFormatter.pad("\(count)", 9) + "  " +1386                    TableFormatter.pad(TableFormatter.fixed(percent, decimals: 2), 9) + "  " +1387                    TableFormatter.pad(TableFormatter.fixed(cumulative, decimals: 2), 9)1388                )1389            }1390            lines.append(String(repeating: "-", count: 15) + "+" + String(repeating: "-", count: 33))1391            lines.append(1392                TableFormatter.pad("Total", 14) + " | " +1393                TableFormatter.pad("\(total)", 9) + "  " +1394                TableFormatter.pad("100.00", 9)1395            )1396            return ZQResult(1397                text: lines.joined(separator: "\n"),1398                scalars: ["N": Double(total), "r": Double(counts.count)]1399            )1400        }14011402        // Two-way table.1403        let second = try labels(names[1])1404        var cells: [String: [String: Int]] = [:]1405        var rowTotals: [String: Int] = [:]1406        var columnTotals: [String: Int] = [:]1407        var total = 01408        for i in 0..<frame.rowCount where mask[i] {1409            let rowLabel = includeMissing ? (first[i] ?? ".") : first[i]1410            let columnLabel = includeMissing ? (second[i] ?? ".") : second[i]1411            guard let rowLabel, let columnLabel else { continue }1412            cells[rowLabel, default: [:]][columnLabel, default: 0] += 11413            rowTotals[rowLabel, default: 0] += 11414            columnTotals[columnLabel, default: 0] += 11415            total += 11416        }1417        guard total > 0 else { throw ZQEngineError("no observations") }14181419        let rows = rowTotals.keys.sorted(by: numericAwareLess)1420        let columns = columnTotals.keys.sorted(by: numericAwareLess)1421        let width = 91422        var lines = [1423            TableFormatter.pad(names[0], 14) + " | " +1424            columns.map { TableFormatter.pad($0, width) }.joined(separator: " ") +1425            "  " + TableFormatter.pad("Total", width),1426            String(repeating: "-", count: 15) + "+" +1427            String(repeating: "-", count: (width + 1) * (columns.count + 1) + 2),1428        ]1429        for row in rows {1430            let cellsText = columns.map {1431                TableFormatter.pad("\(cells[row]?[$0] ?? 0)", width)1432            }.joined(separator: " ")1433            lines.append(1434                TableFormatter.pad(row, 14) + " | " + cellsText + "  " +1435                TableFormatter.pad("\(rowTotals[row] ?? 0)", width)1436            )1437        }1438        lines.append(1439            TableFormatter.pad("Total", 14) + " | " +1440            columns.map { TableFormatter.pad("\(columnTotals[$0] ?? 0)", width) }1441                .joined(separator: " ") +1442            "  " + TableFormatter.pad("\(total)", width)1443        )1444        return ZQResult(1445            text: lines.joined(separator: "\n"),1446            scalars: [1447                "N": Double(total),1448                "r": Double(rows.count),1449                "c": Double(columns.count),1450            ]1451        )1452    }14531454    /// Sorts numeric labels numerically, everything else lexically.1455    private func numericAwareLess(_ a: String, _ b: String) -> Bool {1456        if let x = Double(a), let y = Double(b) { return x < y }1457        return a < b1458    }14591460    private func handleCorrelate(_ command: ZQCommand) throws -> ZQResult {1461        guard !frame.isEmpty else { throw ZQEngineError("no data in memory") }1462        var names = command.varlist.flatMap(\.referencedNames)1463        if names.isEmpty {1464            names = frame.columns.filter(\.data.isNumeric).map(\.name)1465        }1466        guard names.count >= 2 else {1467            throw ZQEngineError("correlate: at least two numeric variables required")1468        }1469        let mask = try observationMask(command)1470        let sources = try names.map { try frame.requireNumeric($0) }14711472        // Listwise deletion across the varlist (Stata correlate).1473        var kept: [Int] = []1474        for i in 0..<frame.rowCount where mask[i] {1475            if sources.allSatisfy({ !$0.missing[i] }) { kept.append(i) }1476        }1477        guard kept.count > 1 else { throw ZQEngineError("no observations") }14781479        let columns = sources.map { source in kept.map { source.values[$0] } }1480        let matrix = ZQCorrelate.matrix(columns: columns)14811482        var lines = [1483            "(obs=\(kept.count))",1484            "",1485            TableFormatter.pad("", 12) + " | " +1486            names.map { TableFormatter.pad($0, 9) }.joined(separator: " "),1487            String(repeating: "-", count: 13) + "+" +1488            String(repeating: "-", count: (9 + 1) * names.count + 1),1489        ]1490        for (i, name) in names.enumerated() {1491            // Lower triangle only, Stata-style.1492            let cells = (0...i).map {1493                TableFormatter.pad(TableFormatter.fixed(matrix[i][$0], decimals: 4), 9)1494            }.joined(separator: " ")1495            lines.append(TableFormatter.pad(name, 12) + " | " + cells)1496        }1497        return ZQResult(1498            text: lines.joined(separator: "\n"),1499            scalars: ["N": Double(kept.count), "rho": matrix[0][1]]1500        )1501    }15021503    private func clusterFromVCE(_ command: ZQCommand) -> String? {1504        guard let vce = command.option("vce"),1505              vce.arguments.count == 2,1506              vce.arguments[0] == "cluster" else { return nil }1507        return vce.arguments[1]1508    }15091510    private func confidenceLevel(_ command: ZQCommand) throws -> Double {1511        guard let text = command.option("level")?.firstArgument else { return 0.95 }1512        guard let value = Double(text), value > 0, value < 100 else {1513            throw ZQEngineError("level() must be between 0 and 100")1514        }1515        return value / 1001516    }15171518    // MARK: - Post-estimation15191520    /// `predict newvar [, xb | residuals | pr | n]` — evaluates the last1521    /// estimation over ALL current observations (missing inputs yield1522    /// missing predictions). Defaults: xb after regress/ivregress, pr1523    /// after logit/probit, n (mean count) after poisson.1524    private func handlePredict(_ command: ZQCommand) throws -> ZQResult {1525        guard let estimation = lastEstimation else {1526            throw ZQEngineError("predict: no estimation results — run a regression first")1527        }1528        guard !frame.isEmpty else { throw ZQEngineError("no data in memory") }1529        let names = command.varlist.flatMap(\.referencedNames)1530        guard names.count == 1, let target = names.first else {1531            throw ZQEngineError("predict: syntax is 'predict newvar [, statistic]'")1532        }15331534        // Boosted models predict through their trees, not a linear xb.1535        if case .boost(let model) = estimation.kind {1536            return try predictBoost(1537                model: model, target: target, command: command,1538                responseName: estimation.responseName1539            )1540        }15411542        // Statistic selection and validation against the model kind.1543        enum Statistic { case xb, residuals, probability, mean }1544        var statistic: Statistic1545        switch estimation.kind {1546        case .ols, .iv: statistic = .xb1547        case .glm(.logit), .glm(.probit): statistic = .probability1548        case .glm(.poisson): statistic = .mean1549        case .boost: fatalError("handled above")1550        }1551        for option in command.options {1552            switch option.name {1553            case "xb": statistic = .xb1554            case "residuals", "resid": statistic = .residuals1555            case "pr":1556                switch estimation.kind {1557                case .glm(.logit), .glm(.probit): statistic = .probability1558                default:1559                    throw ZQEngineError("predict: 'pr' requires logit or probit results")1560                }1561            case "n":1562                guard case .glm(.poisson) = estimation.kind else {1563                    throw ZQEngineError("predict: 'n' requires poisson results")1564                }1565                statistic = .mean1566            default:1567                throw ZQEngineError("predict: unknown statistic '\(option.name)'")1568            }1569        }15701571        // Linear predictor over the full dataset.1572        let n = frame.rowCount1573        var xb = [Double](repeating: 0, count: n)1574        var missing = [Bool](repeating: false, count: n)1575        for coefficient in estimation.coefficients {1576            switch coefficient.definition {1577            case .constant:1578                for i in 0..<n { xb[i] += coefficient.value }1579            case .column(let name):1580                let (values, columnMissing) = try frame.requireNumeric(name)1581                for i in 0..<n {1582                    if columnMissing[i] { missing[i] = true } else {1583                        xb[i] += coefficient.value * values[i]1584                    }1585                }1586            case .indicator(let variable, let level):1587                let (values, columnMissing) = try frame.requireNumeric(variable)1588                for i in 0..<n {1589                    if columnMissing[i] { missing[i] = true }1590                    else if values[i] == level { xb[i] += coefficient.value }1591                }1592            case .product(let variables):1593                let columns = try variables.map { try frame.requireNumeric($0) }1594                for i in 0..<n {1595                    var product = 1.01596                    for column in columns {1597                        if column.missing[i] { missing[i] = true; break }1598                        product *= column.values[i]1599                    }1600                    if !missing[i] { xb[i] += coefficient.value * product }1601                }1602            }1603        }16041605        // Transform.1606        var result = [Double](repeating: .nan, count: n)1607        var label: String1608        switch statistic {1609        case .xb:1610            label = "linear prediction"1611            for i in 0..<n where !missing[i] { result[i] = xb[i] }1612        case .residuals:1613            label = "residuals"1614            let (yValues, yMissing) = try frame.requireNumeric(estimation.responseName)1615            for i in 0..<n where !missing[i] && !yMissing[i] {1616                switch estimation.kind {1617                case .ols, .iv:1618                    result[i] = yValues[i] - xb[i]1619                case .glm(.logit):1620                    result[i] = yValues[i] - 1 / (1 + Foundation.exp(-xb[i]))1621                case .glm(.probit):1622                    result[i] = yValues[i] - ZQDistributions.normalCDF(xb[i])1623                case .glm(.poisson):1624                    result[i] = yValues[i] - Foundation.exp(xb[i])1625                case .boost:1626                    break   // handled by predictBoost1627                }1628            }1629            for i in 0..<n where yMissing[i] { missing[i] = true }1630        case .probability:1631            label = "predicted probability"1632            for i in 0..<n where !missing[i] {1633                if case .glm(.probit) = estimation.kind {1634                    result[i] = ZQDistributions.normalCDF(xb[i])1635                } else {1636                    result[i] = 1 / (1 + Foundation.exp(-xb[i]))1637                }1638            }1639        case .mean:1640            label = "predicted mean count"1641            for i in 0..<n where !missing[i] { result[i] = Foundation.exp(xb[i]) }1642        }1643        for i in 0..<n where missing[i] { result[i] = .nan }16441645        try frame.addColumn(ZQColumn(1646            name: target,1647            data: .float64(values: result, missing: missing)1648        ))1649        let missingCount = missing.count { $0 }1650        var note = "(\(label) '\(target)' generated"1651        if missingCount > 0 { note += ", \(missingCount) missing values" }1652        note += ")"1653        return ZQResult(text: note, scalars: ["N": Double(n - missingCount)])1654    }16551656    /// `elasticnet y x…, lambda(#) [alpha(#)]` and `lasso y x…, lambda(#)`1657    /// — penalized linear regression, glmnet conventions. Without1658    /// lambda(), reports λ_max (the smallest all-zero penalty) as a hint.1659    private func handleElasticNet(1660        _ command: ZQCommand, defaultAlpha: Double?1661    ) throws -> ZQResult {1662        let sample = try buildRegressionSample(command, clusterVariable: nil)16631664        var alpha = defaultAlpha ?? 11665        if let text = command.option("alpha")?.firstArgument {1666            guard defaultAlpha == nil else {1667                throw ZQEngineError("lasso: alpha is fixed at 1 — use elasticnet for other mixes")1668            }1669            guard let value = Double(text) else {1670                throw ZQEngineError("elasticnet: invalid alpha()")1671            }1672            alpha = value1673        }1674        guard let lambdaText = command.option("lambda")?.firstArgument,1675              let lambda = Double(lambdaText) else {1676            let hint = ZQElasticNet.lambdaMax(1677                y: sample.y, predictors: sample.predictors, alpha: alpha1678            )1679            throw ZQEngineError(1680                "\(command.verb): lambda(#) required — lambda_max for this data is \(TableFormatter.general(hint, significant: 6))"1681            )1682        }16831684        let result = try ZQElasticNet.fit(1685            y: sample.y,1686            predictors: sample.predictors,1687            alpha: alpha,1688            lambda: lambda1689        )16901691        // predict works after penalized fits; margins (needing a VCE)1692        // deliberately does not.1693        recordEstimation(1694            kind: .ols,1695            responseName: sample.responseName,1696            coefficients: result.coefficients.map {1697                ZQCoefficient(1698                    name: $0.name, estimate: $0.value, standardError: .nan,1699                    tStatistic: .nan, pValue: .nan,1700                    confidenceLower: .nan, confidenceUpper: .nan1701                )1702            } + [ZQCoefficient(1703                name: "_cons", estimate: result.intercept, standardError: .nan,1704                tStatistic: .nan, pValue: .nan,1705                confidenceLower: .nan, confidenceUpper: .nan1706            )],1707            definitions: sample.definitions,1708            includeConstant: true,1709            vce: [],1710            sampleSize: sample.y.count1711        )17121713        var lines = [1714            alpha == 1 ? "Lasso linear model" : "Elastic-net linear model (alpha = \(TableFormatter.general(alpha)))",1715        ]1716        if sample.droppedMissing > 0 {1717            lines.append("(\(sample.droppedMissing) observations dropped due to missing values)")1718        }1719        for (label, value) in [1720            ("Number of obs", "\(sample.y.count)"),1721            ("lambda", TableFormatter.general(lambda, significant: 6)),1722            ("Nonzero coefficients", "\(result.nonzeroCount) of \(result.coefficients.count)"),1723            ("R-squared", TableFormatter.fixed(result.rSquared, decimals: 4)),1724        ] {1725            lines.append(1726                TableFormatter.pad(label, 46, right: false) + "= " +1727                TableFormatter.pad(value, 12)1728            )1729        }1730        lines.append("")1731        lines.append(1732            TableFormatter.pad(sample.responseName, 12) + " | " +1733            TableFormatter.pad("Coefficient", 12)1734        )1735        lines.append(String(repeating: "-", count: 13) + "+" + String(repeating: "-", count: 14))17361737        var scalars: [String: Double] = [1738            "N": Double(sample.y.count),1739            "lambda": lambda,1740            "alpha": alpha,1741            "k_nonzero": Double(result.nonzeroCount),1742            "r2": result.rSquared,1743        ]1744        for coefficient in result.coefficients where coefficient.value != 0 {1745            lines.append(1746                TableFormatter.pad(coefficient.name, 12) + " | " +1747                TableFormatter.pad(TableFormatter.general(coefficient.value), 12)1748            )1749            scalars["b_\(coefficient.name)"] = coefficient.value1750        }1751        // Zeroed coefficients still surface as scalars (as exact 0).1752        for coefficient in result.coefficients where coefficient.value == 0 {1753            scalars["b_\(coefficient.name)"] = 01754        }1755        lines.append(1756            TableFormatter.pad("_cons", 12) + " | " +1757            TableFormatter.pad(TableFormatter.general(result.intercept), 12)1758        )1759        scalars["b__cons"] = result.intercept1760        return ZQResult(text: lines.joined(separator: "\n"), scalars: scalars)1761    }17621763    /// `boost y x…, rounds(#) [eta(#) maxdepth(#) lambda(#)]` —1764    /// gradient-boosted regression trees (xgboost exact-greedy, squared1765    /// loss). Deterministic: no subsampling.1766    private func handleBoost(_ command: ZQCommand) throws -> ZQResult {1767        let sample = try buildRegressionSample(command, clusterVariable: nil)1768        func numberOption(_ name: String, default defaultValue: Double) throws -> Double {1769            guard let text = command.option(name)?.firstArgument else { return defaultValue }1770            guard let value = Double(text), value > 0 else {1771                throw ZQEngineError("boost: invalid \(name)()")1772            }1773            return value1774        }1775        guard let roundsText = command.option("rounds")?.firstArgument,1776              let rounds = Int(roundsText), rounds > 0 else {1777            throw ZQEngineError("boost: rounds(#) required")1778        }1779        let model = try ZQGradientBoosting.fit(1780            y: sample.y,1781            features: sample.predictors,1782            rounds: rounds,1783            learningRate: try numberOption("eta", default: 0.3),1784            maxDepth: Int(try numberOption("maxdepth", default: 6)),1785            lambda: try numberOption("lambda", default: 1)1786        )17871788        // Training fit.1789        let n = sample.y.count1790        var rss = 0.01791        var tss = 0.01792        let meanY = sample.y.reduce(0, +) / Double(n)1793        for i in 0..<n {1794            let vector = sample.predictors.map { $0.values[i] }1795            let predicted = model.predict(vector)1796            rss += (sample.y[i] - predicted) * (sample.y[i] - predicted)1797            tss += (sample.y[i] - meanY) * (sample.y[i] - meanY)1798        }1799        let rSquared = tss > 0 ? 1 - rss / tss : .nan1800        let rmse = (rss / Double(n)).squareRoot()18011802        recordEstimation(1803            kind: .boost(model),1804            responseName: sample.responseName,1805            coefficients: sample.predictors.map {1806                ZQCoefficient(1807                    name: $0.name, estimate: .nan, standardError: .nan,1808                    tStatistic: .nan, pValue: .nan,1809                    confidenceLower: .nan, confidenceUpper: .nan1810                )1811            },1812            definitions: sample.definitions,1813            includeConstant: false,1814            sampleSize: n1815        )18161817        var lines = ["Gradient-boosted trees (squared loss)"]1818        if sample.droppedMissing > 0 {1819            lines.append("(\(sample.droppedMissing) observations dropped due to missing values)")1820        }1821        for (label, value) in [1822            ("Number of obs", "\(n)"),1823            ("Trees", "\(rounds)"),1824            ("Learning rate (eta)", TableFormatter.general(model.learningRate)),1825            ("Max depth", command.option("maxdepth")?.firstArgument ?? "6"),1826            ("Training R-squared", TableFormatter.fixed(rSquared, decimals: 4)),1827            ("Training RMSE", TableFormatter.general(rmse, significant: 6)),1828        ] {1829            lines.append(1830                TableFormatter.pad(label, 46, right: false) + "= " +1831                TableFormatter.pad(value, 12)1832            )1833        }1834        lines.append("")1835        lines.append("Use 'predict newvar' for fitted values. Training fit is in-sample —")1836        lines.append("expect it to be optimistic relative to held-out data.")1837        return ZQResult(1838            text: lines.joined(separator: "\n"),1839            scalars: [1840                "N": Double(n), "rounds": Double(rounds),1841                "r2": rSquared, "rmse": rmse,1842            ]1843        )1844    }18451846    /// `margins, dydx(varlist)` — average marginal effects with1847    /// delta-method standard errors. OLS/IV effects are the coefficients1848    /// themselves; GLM effects average dμ/dx over the estimation sample:1849    ///   AME_j = β_j · (1/N) Σᵢ g′(xbᵢ)1850    ///   ∂AME_j/∂β_m = (1/N) Σᵢ g″(xbᵢ)·x_im·β_j + δ_jm·(1/N) Σᵢ g′(xbᵢ)1851    /// Only continuous (plain-column) regressors are supported so far.1852    private func handleMargins(_ command: ZQCommand) throws -> ZQResult {1853        guard let estimation = lastEstimation else {1854            throw ZQEngineError("margins: no estimation results — run a regression first")1855        }1856        guard let dydx = command.option("dydx"), !dydx.arguments.isEmpty else {1857            throw ZQEngineError("margins: syntax is 'margins, dydx(varlist)'")1858        }1859        guard !estimation.vce.isEmpty else {1860            throw ZQEngineError("margins: the last estimation stored no covariance matrix")1861        }18621863        let k = estimation.coefficients.count1864        struct Effect {1865            var name: String1866            var dydx: Double1867            var standardError: Double1868        }1869        var effects: [Effect] = []18701871        for name in dydx.arguments {1872            guard let j = estimation.coefficients.firstIndex(where: {1873                $0.definition == .column(name)1874            }) else {1875                if estimation.coefficients.contains(where: {1876                    if case .indicator(let variable, _) = $0.definition { return variable == name }1877                    if case .product(let variables) = $0.definition { return variables.contains(name) }1878                    return false1879                }) {1880                    throw ZQEngineError(1881                        "margins: dydx over factor or interaction terms is not supported yet"1882                    )1883                }1884                throw ZQEngineError("margins: '\(name)' is not a regressor in the last estimation")1885            }1886            let beta = estimation.coefficients[j].value18871888            switch estimation.kind {1889            case .boost:1890                throw ZQEngineError("margins: not available after boost")1891            case .ols, .iv:1892                // Linear model: the marginal effect is the coefficient.1893                effects.append(Effect(1894                    name: name,1895                    dydx: beta,1896                    standardError: estimation.vce[j * k + j].squareRoot()1897                ))18981899            case .glm(let family):1900                guard let design = estimation.sampleDesign,1901                      estimation.sampleSize > 0 else {1902                    throw ZQEngineError("margins: estimation sample unavailable")1903                }1904                let n = estimation.sampleSize1905                // xb over the estimation sample.1906                var xb = [Double](repeating: 0, count: n)1907                for m in 0..<k {1908                    let value = estimation.coefficients[m].value1909                    for i in 0..<n { xb[i] += value * design[m * n + i] }1910                }1911                // g′ and g″ per observation.1912                var meanFirst = 0.01913                var second = [Double](repeating: 0, count: n)1914                for i in 0..<n {1915                    switch family {1916                    case .logit:1917                        let p = 1 / (1 + Foundation.exp(-xb[i]))1918                        meanFirst += p * (1 - p)1919                        second[i] = p * (1 - p) * (1 - 2 * p)1920                    case .probit:1921                        let density = Foundation.exp(-0.5 * xb[i] * xb[i])1922                            / (2 * Double.pi).squareRoot()1923                        meanFirst += density1924                        second[i] = -xb[i] * density1925                    case .poisson:1926                        let mu = Foundation.exp(xb[i])1927                        meanFirst += mu1928                        second[i] = mu1929                    }1930                }1931                meanFirst /= Double(n)19321933                // Delta-method gradient.1934                var gradient = [Double](repeating: 0, count: k)1935                for m in 0..<k {1936                    var meanSecondX = 0.01937                    for i in 0..<n { meanSecondX += second[i] * design[m * n + i] }1938                    gradient[m] = beta * meanSecondX / Double(n)1939                }1940                gradient[j] += meanFirst19411942                var variance = 0.01943                for a in 0..<k {1944                    for b in 0..<k {1945                        variance += gradient[a] * estimation.vce[b * k + a] * gradient[b]1946                    }1947                }1948                effects.append(Effect(1949                    name: name,1950                    dydx: beta * meanFirst,1951                    standardError: variance.squareRoot()1952                ))1953            }1954        }19551956        // Render: t inference when the model carries a df, z otherwise.1957        let level = try confidenceLevel(command)1958        let usesT = estimation.inferenceDF != nil1959        let critical: Double1960        if let df = estimation.inferenceDF {1961            critical = ZQDistributions.studentTQuantile(0.5 + level / 2, df: df)1962        } else {1963            critical = ZQDistributions.normalQuantile(0.5 + level / 2)1964        }19651966        var lines = [1967            "Average marginal effects                      Number of obs = " +1968            TableFormatter.pad("\(estimation.sampleSize)", 10),1969            "",1970        ]1971        let widths = [12, 12, 11, 8, 8, 22]1972        let statLabel = usesT ? "t" : "z"1973        lines.append(1974            TableFormatter.pad("", widths[0]) + " | " +1975            TableFormatter.pad("dy/dx", widths[1]) + "  " +1976            TableFormatter.pad("Std. err.", widths[2]) + "  " +1977            TableFormatter.pad(statLabel, widths[3]) + "  " +1978            TableFormatter.pad("P>|\(statLabel)|", widths[4]) + "  " +1979            TableFormatter.pad("[\(Int(level * 100))% conf. interval]", widths[5])1980        )1981        lines.append(TableFormatter.rule(widths))19821983        var scalars: [String: Double] = ["N": Double(estimation.sampleSize)]1984        for effect in effects {1985            let statistic = effect.dydx / effect.standardError1986            let p: Double1987            if let df = estimation.inferenceDF {1988                p = ZQDistributions.tTestPValue(statistic, df: df)1989            } else {1990                p = 2 * (1 - ZQDistributions.normalCDF(abs(statistic)))1991            }1992            lines.append(1993                TableFormatter.pad(effect.name, widths[0]) + " | " +1994                TableFormatter.pad(TableFormatter.general(effect.dydx), widths[1]) + "  " +1995                TableFormatter.pad(TableFormatter.general(effect.standardError), widths[2]) + "  " +1996                TableFormatter.pad(TableFormatter.fixed(statistic, decimals: 2), widths[3]) + "  " +1997                TableFormatter.pad(TableFormatter.fixed(p, decimals: 3), widths[4]) + "  " +1998                TableFormatter.pad(1999                    TableFormatter.general(effect.dydx - critical * effect.standardError) + "   " +2000                    TableFormatter.general(effect.dydx + critical * effect.standardError),2001                    widths[5]2002                )2003            )2004            scalars["dydx_\(effect.name)"] = effect.dydx2005            scalars["se_\(effect.name)"] = effect.standardError2006        }2007        return ZQResult(text: lines.joined(separator: "\n"), scalars: scalars)2008    }20092010    /// Boosted-tree predictions: trees route missing features down the2011    /// default (left) branch, so predictions exist for every observation.2012    private func predictBoost(2013        model: ZQBoostModel, target: String, command: ZQCommand,2014        responseName: String2015    ) throws -> ZQResult {2016        for option in command.options2017        where !["xb", "residuals", "resid"].contains(option.name) {2018            throw ZQEngineError("predict: '\(option.name)' is not available after boost")2019        }2020        let wantsResiduals = command.hasOption("residuals") || command.hasOption("resid")20212022        let n = frame.rowCount2023        let columns = try model.featureNames.map { try frame.requireNumeric($0) }2024        var values = [Double](repeating: .nan, count: n)2025        var missing = [Bool](repeating: false, count: n)2026        var yValues = [Double](repeating: 0, count: n)2027        var yMissing = [Bool](repeating: false, count: n)2028        if wantsResiduals {2029            (yValues, yMissing) = try frame.requireNumeric(responseName)2030        }20312032        for i in 0..<n {2033            let vector = columns.map { $0.missing[i] ? Double.nan : $0.values[i] }2034            let predicted = model.predict(vector)2035            if wantsResiduals {2036                if yMissing[i] {2037                    missing[i] = true2038                } else {2039                    values[i] = yValues[i] - predicted2040                }2041            } else {2042                values[i] = predicted2043            }2044        }2045        try frame.addColumn(ZQColumn(2046            name: target, data: .float64(values: values, missing: missing)2047        ))2048        let label = wantsResiduals ? "residuals" : "boosted prediction"2049        return ZQResult(text: "(\(label) '\(target)' generated)")2050    }20512052    // MARK: - Bootstrap prefix20532054    private func handleBootstrap(2055        _ command: ZQCommand, backend: ZQBackend2056    ) async throws -> ZQResult {2057        guard let body = command.body else {2058            throw ZQEngineError("bootstrap: syntax is 'bootstrap, reps(#): command'")2059        }2060        guard body.verb == "regress" else {2061            throw ZQEngineError("bootstrap currently supports 'regress' bodies only")2062        }2063        guard let repsText = command.option("reps")?.firstArgument,2064              let reps = Int(repsText), reps >= 2 else {2065            throw ZQEngineError("bootstrap: reps(#) with # ≥ 2 required")2066        }2067        if let seedText = command.option("seed")?.firstArgument {2068            guard let value = UInt64(seedText) else {2069                throw ZQEngineError("bootstrap: invalid seed")2070            }2071            seed = value2072        }20732074        let sample = try buildRegressionSample(body, clusterVariable: nil)2075        let includeConstant = !body.hasOption("noconstant")2076        let point = try ZQOLS.fit(2077            y: sample.y,2078            predictors: sample.predictors,2079            includeConstant: includeConstant,2080            variance: .classical2081        )20822083        let n = sample.y.count2084        let k = point.coefficients.count2085        let generator = Philox4x32(seed: seed)2086        let y = sample.y2087        let predictors = sample.predictors20882089        // Pairs bootstrap: replicate r regenerates its indices from the2090        // Philox counter stream — deterministic and order-independent, so2091        // replicates can run in parallel CPU chunks or as GPU batches with2092        // identical resamples.2093        let allDraws: [[Double]]2094        if backend == .gpu {2095            allDraws = ZQGPUBootstrap.pairsBootstrapOLS(2096                y: y,2097                predictors: predictors,2098                includeConstant: includeConstant,2099                replicates: reps,2100                seed: seed2101            )2102        } else {2103            allDraws = try await withThrowingTaskGroup(2104                of: [(Int, [Double])].self2105            ) { group in2106                let chunkCount = min(2107                    max(1, ProcessInfo.processInfo.activeProcessorCount), reps2108                )2109                let chunkSize = (reps + chunkCount - 1) / chunkCount2110                for chunk in 0..<chunkCount {2111                    let lower = chunk * chunkSize2112                    let upper = min(reps, lower + chunkSize)2113                    guard lower < upper else { continue }2114                    group.addTask {2115                        var results: [(Int, [Double])] = []2116                        results.reserveCapacity(upper - lower)2117                        for replicate in lower..<upper {2118                            let indices = ZQResampling.pairsBootstrapIndices(2119                                replicate: replicate, sampleSize: n, generator: generator2120                            )2121                            let resampledY = indices.map { y[$0] }2122                            let resampledPredictors = predictors.map { column in2123                                (name: column.name, values: indices.map { column.values[$0] })2124                            }2125                            let fit = try ZQOLS.fit(2126                                y: resampledY,2127                                predictors: resampledPredictors,2128                                includeConstant: includeConstant,2129                                variance: .classical2130                            )2131                            results.append((replicate, fit.coefficients.map(\.estimate)))2132                        }2133                        return results2134                    }2135                }2136                var collected = [[Double]](repeating: [], count: reps)2137                for try await chunkResults in group {2138                    for (replicate, estimates) in chunkResults {2139                        collected[replicate] = estimates2140                    }2141                }2142                return collected2143            }2144        }21452146        // Singular resamples (NaN rows from the GPU path) are dropped,2147        // Stata-style.2148        let draws = allDraws.filter { row in !row.contains(where: \.isNaN) }2149        let completed = draws.count2150        guard completed >= 2 else {2151            throw ZQEngineError("bootstrap: too many failed replicates")2152        }21532154        // Bootstrap standard errors: sd of replicate estimates.2155        let title = backend == .gpu2156            ? "Bootstrap results (pairs, GPU batched)"2157            : "Bootstrap results (pairs)"2158        var lines = [2159            TableFormatter.pad(title, 46, right: false) + "Replications  = " +2160            TableFormatter.pad("\(completed)", 10),2161            "                                              Number of obs = " +2162            TableFormatter.pad("\(n)", 10),2163        ]2164        if completed < reps {2165            lines.append("(\(reps - completed) replicates dropped: singular resample)")2166        }2167        lines.append("")2168        let widths = [12, 12, 12, 8, 8, 22]2169        lines.append(2170            TableFormatter.pad(sample.responseName, widths[0]) + " | " +2171            TableFormatter.pad("Observed", widths[1]) + "  " +2172            TableFormatter.pad("Bootstrap", widths[2]) + "  " +2173            TableFormatter.pad("z", widths[3]) + "  " +2174            TableFormatter.pad("P>|z|", widths[4]) + "  " +2175            TableFormatter.pad("[95% conf. interval]", widths[5])2176        )2177        lines.append(2178            TableFormatter.pad("", widths[0]) + " | " +2179            TableFormatter.pad("coefficient", widths[1]) + "  " +2180            TableFormatter.pad("std. err.", widths[2])2181        )2182        lines.append(TableFormatter.rule(widths))21832184        var scalars: [String: Double] = [2185            "reps": Double(completed),2186            "N": Double(n),2187            "gpu": backend == .gpu ? 1 : 0,2188        ]2189        for j in 0..<k {2190            let estimates = draws.map { $0[j] }2191            let mean = estimates.reduce(0, +) / Double(completed)2192            let variance = estimates.reduce(0) { $0 + ($1 - mean) * ($1 - mean) }2193                / Double(completed - 1)2194            let se = variance.squareRoot()2195            let observed = point.coefficients[j].estimate2196            let z = observed / se2197            let p = 2 * (1 - ZQDistributions.normalCDF(abs(z)))2198            let critical = 1.959963984540054  // Φ⁻¹(0.975)2199            lines.append(2200                TableFormatter.pad(point.coefficients[j].name, widths[0]) + " | " +2201                TableFormatter.pad(TableFormatter.general(observed), widths[1]) + "  " +2202                TableFormatter.pad(TableFormatter.general(se), widths[2]) + "  " +2203                TableFormatter.pad(TableFormatter.fixed(z, decimals: 2), widths[3]) + "  " +2204                TableFormatter.pad(TableFormatter.fixed(p, decimals: 3), widths[4]) + "  " +2205                TableFormatter.pad(2206                    TableFormatter.general(observed - critical * se) + "   " +2207                    TableFormatter.general(observed + critical * se),2208                    widths[5]2209                )2210            )2211            scalars["b_\(point.coefficients[j].name)"] = observed2212            scalars["se_\(point.coefficients[j].name)"] = se2213        }2214        return ZQResult(text: lines.joined(separator: "\n"), scalars: scalars)2215    }22162217    // MARK: - Bayesian estimation22182219    /// `bayes [, mcmcsize(#) burnin(#) seed(#) normalprior(#)]: reg y x…`2220    /// — Gibbs-sampled Bayesian linear regression with Stata-style default2221    /// priors: N(0, 10000) on coefficients, InvGamma(0.01, 0.01) on σ².2222    private func handleBayes(_ command: ZQCommand) throws -> ZQResult {2223        guard let body = command.body else {2224            throw ZQEngineError("bayes: syntax is 'bayes [, options]: regress …'")2225        }2226        guard body.verb == "regress" else {2227            throw ZQEngineError("bayes currently supports 'regress' bodies only")2228        }2229        if let seedText = command.option("seed")?.firstArgument {2230            guard let value = UInt64(seedText) else {2231                throw ZQEngineError("bayes: invalid seed")2232            }2233            seed = value2234        }2235        func intOption(_ name: String, default defaultValue: Int) throws -> Int {2236            guard let text = command.option(name)?.firstArgument else { return defaultValue }2237            guard let value = Int(text), value > 0 else {2238                throw ZQEngineError("bayes: invalid \(name)()")2239            }2240            return value2241        }2242        let mcmcSize = try intOption("mcmcsize", default: 10_000)2243        let burnIn = try intOption("burnin", default: 2_500)2244        var priorVariance = 10_000.02245        if let text = command.option("normalprior")?.firstArgument {2246            guard let value = Double(text), value > 0 else {2247                throw ZQEngineError("bayes: invalid normalprior()")2248            }2249            priorVariance = value2250        }22512252        let sample = try buildRegressionSample(body, clusterVariable: nil)2253        let result = try ZQBayesianRegression.fitGibbs(2254            y: sample.y,2255            predictors: sample.predictors,2256            includeConstant: !body.hasOption("noconstant"),2257            mcmcSize: mcmcSize,2258            burnIn: burnIn,2259            seed: seed,2260            coefficientPriorVariance: priorVariance2261        )2262        lastEstimation = nil   // predict after bayes needs posterior draws22632264        var lines = ["Bayesian linear regression (Gibbs)"]2265        if sample.droppedMissing > 0 {2266            lines.append("(\(sample.droppedMissing) observations dropped due to missing values)")2267        }2268        for (label, value) in [2269            ("Number of obs", "\(result.observationCount)"),2270            ("MCMC iterations", "\(result.mcmcSize)"),2271            ("Burn-in", "\(result.burnIn)"),2272            ("Priors", "b ~ N(0, \(TableFormatter.general(priorVariance))), sigma2 ~ IG(.01, .01)"),2273        ] {2274            lines.append(2275                TableFormatter.pad(label, 46, right: false) + "= " +2276                TableFormatter.pad(value, 24)2277            )2278        }2279        lines.append("")22802281        let widths = [12, 12, 12, 24]2282        lines.append(2283            TableFormatter.pad(sample.responseName, widths[0]) + " | " +2284            TableFormatter.pad("Mean", widths[1]) + "  " +2285            TableFormatter.pad("Std. dev.", widths[2]) + "  " +2286            TableFormatter.pad("[95% cred. interval]", widths[3])2287        )2288        lines.append(TableFormatter.rule(widths))22892290        var scalars: [String: Double] = [2291            "N": Double(result.observationCount),2292            "mcmcsize": Double(result.mcmcSize),2293            "burnin": Double(result.burnIn),2294        ]2295        for coefficient in result.coefficients + [result.sigma] {2296            lines.append(2297                TableFormatter.pad(coefficient.name, widths[0]) + " | " +2298                TableFormatter.pad(TableFormatter.general(coefficient.posteriorMean), widths[1]) + "  " +2299                TableFormatter.pad(TableFormatter.general(coefficient.posteriorSD), widths[2]) + "  " +2300                TableFormatter.pad(2301                    TableFormatter.general(coefficient.credibleLower) + "   " +2302                    TableFormatter.general(coefficient.credibleUpper),2303                    widths[3]2304                )2305            )2306            scalars["b_\(coefficient.name)"] = coefficient.posteriorMean2307            scalars["sd_\(coefficient.name)"] = coefficient.posteriorSD2308        }2309        return ZQResult(text: lines.joined(separator: "\n"), scalars: scalars)2310    }23112312    // MARK: - Permutation test23132314    /// `permute, reps(#) [seed(#)]: reg y x…` — permutes the response2315    /// within the estimation sample, refits, and reports the empirical2316    /// two-sided p-value per coefficient: the fraction of permuted |β*|2317    /// at or above the observed |β|. Permutations come from the Philox2318    /// argsort stream, so results are exactly reproducible and2319    /// replicate-addressable.2320    private func handlePermute(_ command: ZQCommand) async throws -> ZQResult {2321        guard let body = command.body else {2322            throw ZQEngineError("permute: syntax is 'permute, reps(#): command'")2323        }2324        guard body.verb == "regress" else {2325            throw ZQEngineError("permute currently supports 'regress' bodies only")2326        }2327        guard let repsText = command.option("reps")?.firstArgument,2328              let reps = Int(repsText), reps >= 1 else {2329            throw ZQEngineError("permute: reps(#) with # ≥ 1 required")2330        }2331        if let seedText = command.option("seed")?.firstArgument {2332            guard let value = UInt64(seedText) else {2333                throw ZQEngineError("permute: invalid seed")2334            }2335            seed = value2336        }23372338        let sample = try buildRegressionSample(body, clusterVariable: nil)2339        let includeConstant = !body.hasOption("noconstant")2340        let observed = try ZQOLS.fit(2341            y: sample.y,2342            predictors: sample.predictors,2343            includeConstant: includeConstant,2344            variance: .classical2345        )2346        let n = sample.y.count2347        let k = observed.coefficients.count2348        let generator = Philox4x32(seed: seed)2349        let y = sample.y2350        let predictors = sample.predictors23512352        // Count |β*| ≥ |β_obs| per coefficient, in parallel chunks —2353        // permutations are counter-addressable so order is irrelevant.2354        let observedMagnitudes = observed.coefficients.map { abs($0.estimate) }2355        let chunkCount = min(2356            max(1, ProcessInfo.processInfo.activeProcessorCount), reps2357        )2358        let chunkSize = (reps + chunkCount - 1) / chunkCount2359        let exceedances: [Int] = try await withThrowingTaskGroup(of: [Int].self) { group in2360            for chunk in 0..<chunkCount {2361                let lower = chunk * chunkSize2362                let upper = min(reps, lower + chunkSize)2363                guard lower < upper else { continue }2364                group.addTask {2365                    var counts = [Int](repeating: 0, count: k)2366                    for replicate in lower..<upper {2367                        let order = ZQResampling.permutationIndices(2368                            replicate: replicate, sampleSize: n, generator: generator2369                        )2370                        let fit = try ZQOLS.fit(2371                            y: order.map { y[$0] },2372                            predictors: predictors,2373                            includeConstant: includeConstant,2374                            variance: .classical2375                        )2376                        for j in 0..<k2377                        where abs(fit.coefficients[j].estimate) >= observedMagnitudes[j] {2378                            counts[j] += 12379                        }2380                    }2381                    return counts2382                }2383            }2384            var total = [Int](repeating: 0, count: k)2385            for try await counts in group {2386                for j in 0..<k { total[j] += counts[j] }2387            }2388            return total2389        }23902391        var lines = [2392            TableFormatter.pad("Permutation test (response permuted)", 46, right: false) +2393            "Replications  = " + TableFormatter.pad("\(reps)", 10),2394            "                                              Number of obs = " +2395            TableFormatter.pad("\(n)", 10),2396            "",2397        ]2398        let widths = [12, 12, 8, 8, 10, 11]2399        lines.append(2400            TableFormatter.pad(sample.responseName, widths[0]) + " | " +2401            TableFormatter.pad("Observed", widths[1]) + "  " +2402            TableFormatter.pad("c", widths[2]) + "  " +2403            TableFormatter.pad("reps", widths[3]) + "  " +2404            TableFormatter.pad("p=c/reps", widths[4]) + "  " +2405            TableFormatter.pad("SE(p)", widths[5])2406        )2407        lines.append(TableFormatter.rule(widths))24082409        var scalars: [String: Double] = ["reps": Double(reps), "N": Double(n)]2410        for j in 0..<k {2411            let name = observed.coefficients[j].name2412            // The constant's permutation distribution is degenerate under2413            // response permutation; Stata omits it too.2414            if name == "_cons" { continue }2415            let p = Double(exceedances[j]) / Double(reps)2416            let standardError = (p * (1 - p) / Double(reps)).squareRoot()2417            lines.append(2418                TableFormatter.pad(name, widths[0]) + " | " +2419                TableFormatter.pad(2420                    TableFormatter.general(observed.coefficients[j].estimate), widths[1]2421                ) + "  " +2422                TableFormatter.pad("\(exceedances[j])", widths[2]) + "  " +2423                TableFormatter.pad("\(reps)", widths[3]) + "  " +2424                TableFormatter.pad(TableFormatter.fixed(p, decimals: 4), widths[4]) + "  " +2425                TableFormatter.pad(TableFormatter.fixed(standardError, decimals: 4), widths[5])2426            )2427            scalars["b_\(name)"] = observed.coefficients[j].estimate2428            scalars["p_\(name)"] = p2429            scalars["c_\(name)"] = Double(exceedances[j])2430        }2431        return ZQResult(text: lines.joined(separator: "\n"), scalars: scalars)2432    }24332434    // MARK: - Graphics24352436    private func handleGraph(_ command: ZQCommand) throws -> ZQResult {2437        guard !frame.isEmpty else { throw ZQEngineError("no data in memory") }2438        let kind: ZQPlotSpec.Kind2439        switch command.subverb {2440        case "scatter": kind = .scatter2441        case "line": kind = .line2442        case "histogram": kind = .histogram2443        case "kdensity": kind = .kdensity2444        case let other:2445            throw ZQEngineError("graph \(other ?? ""): not supported yet")2446        }24472448        let names = command.varlist.flatMap(\.referencedNames)2449        let mask = try observationMask(command)24502451        // Single-variable plots: histogram and kernel density.2452        if kind == .histogram || kind == .kdensity {2453            guard names.count == 1 else {2454                throw ZQEngineError("\(command.subverb ?? ""): syntax is '\(command.subverb ?? "") varname'")2455            }2456            let (values, missing) = try frame.requireNumeric(names[0])2457            let data = (0..<frame.rowCount)2458                .filter { mask[$0] && !missing[$0] }2459                .map { values[$0] }2460            guard data.count > 1 else { throw ZQEngineError("no observations") }24612462            let series: ZQPlotSpec.Series2463            if kind == .histogram {2464                series = Self.histogramSeries(2465                    data: data, label: names[0],2466                    binCount: command.option("bins")?.firstArgument.flatMap(Int.init)2467                )2468            } else {2469                series = Self.kernelDensitySeries(data: data, label: names[0])2470            }2471            lastPlot = ZQPlotSpec(2472                kind: kind,2473                xLabel: names[0],2474                yLabel: kind == .histogram ? "Frequency" : "Density",2475                series: [series]2476            )2477            return ZQResult(2478                text: "(plot created: \(data.count) observations)",2479                scalars: ["N": Double(data.count)]2480            )2481        }24822483        guard names.count == 2 else {2484            throw ZQEngineError("graph \(command.subverb ?? ""): syntax is 'graph \(command.subverb ?? "kind") yvar xvar'")2485        }2486        let (yValues, yMissing) = try frame.requireNumeric(names[0])2487        let (xValues, xMissing) = try frame.requireNumeric(names[1])24882489        func series(label: String, rows: [Int]) -> ZQPlotSpec.Series {2490            ZQPlotSpec.Series(2491                label: label,2492                x: rows.map { xValues[$0] },2493                y: rows.map { yValues[$0] }2494            )2495        }24962497        let complete = (0..<frame.rowCount).filter {2498            mask[$0] && !yMissing[$0] && !xMissing[$0]2499        }25002501        var allSeries: [ZQPlotSpec.Series] = []2502        if let byName = command.option("by")?.firstArgument {2503            let (byValues, byMissing) = try frame.requireNumeric(byName)2504            let levels = ZQFactorExpansion.levels(2505                values: complete.map { byValues[$0] },2506                missing: complete.map { byMissing[$0] }2507            )2508            for level in levels {2509                let rows = complete.filter { !byMissing[$0] && byValues[$0] == level }2510                let label = level == level.rounded()2511                    ? "\(byName)=\(Int(level))"2512                    : "\(byName)=\(level)"2513                allSeries.append(series(label: label, rows: rows))2514            }2515        } else {2516            allSeries.append(series(label: names[0], rows: complete))2517        }25182519        lastPlot = ZQPlotSpec(2520            kind: kind,2521            xLabel: names[1],2522            yLabel: names[0],2523            series: allSeries2524        )2525        let pointCount = allSeries.reduce(0) { $0 + $1.x.count }2526        return ZQResult(text: "(plot created: \(pointCount) points)")2527    }25282529    // MARK: - UI support25302531    /// Compiles filter-bar text into a keep-mask over the current dataset2532    /// (the data browser's filter compiles to an `if` expression, §7).2533    /// Missing conditions exclude the observation, like `if` qualifiers.2534    public func conditionMask(_ expressionText: String) throws -> [Bool] {2535        let expression = try parser.parseExpression(expressionText)2536        let evaluator = ExpressionEvaluator(frame: frame)2537        return try evaluator.evaluateCondition(expression)2538    }25392540    /// Frequency histogram: Sturges bin count by default, half-open bins2541    /// with the maximum folded into the last bin. Series points are bin2542    /// midpoints vs counts.2543    static func histogramSeries(2544        data: [Double], label: String, binCount: Int?2545    ) -> ZQPlotSpec.Series {2546        let n = data.count2547        let bins = max(1, binCount ?? Int((Foundation.log2(Double(n))).rounded(.up)) + 1)2548        let minimum = data.min()!2549        let maximum = data.max()!2550        let width = maximum > minimum ? (maximum - minimum) / Double(bins) : 12551        var counts = [Double](repeating: 0, count: bins)2552        for value in data {2553            let raw = Int((value - minimum) / width)2554            counts[min(max(raw, 0), bins - 1)] += 12555        }2556        let midpoints = (0..<bins).map { minimum + (Double($0) + 0.5) * width }2557        return ZQPlotSpec.Series(label: label, x: midpoints, y: counts)2558    }25592560    /// Epanechnikov kernel density on a 100-point grid with Silverman's2561    /// bandwidth h = 0.9·min(sd, IQR/1.349)·n^(−1/5) (Stata's default).2562    static func kernelDensitySeries(2563        data: [Double], label: String2564    ) -> ZQPlotSpec.Series {2565        let n = data.count2566        let mean = data.reduce(0, +) / Double(n)2567        let sd = (data.reduce(0) { $0 + ($1 - mean) * ($1 - mean) }2568            / Double(n - 1)).squareRoot()2569        let sorted = data.sorted()2570        let iqr = ZQSummarize.stataPercentile(sorted: sorted, percent: 75)2571            - ZQSummarize.stataPercentile(sorted: sorted, percent: 25)2572        var spread = min(sd, iqr / 1.349)2573        if spread <= 0 { spread = max(sd, 1e-9) }2574        let bandwidth = 0.9 * spread * pow(Double(n), -0.2)25752576        let lower = sorted.first! - bandwidth2577        let upper = sorted.last! + bandwidth2578        let gridSize = 1002579        let step = (upper - lower) / Double(gridSize - 1)2580        var xs = [Double](repeating: 0, count: gridSize)2581        var densities = [Double](repeating: 0, count: gridSize)2582        for g in 0..<gridSize {2583            let point = lower + Double(g) * step2584            xs[g] = point2585            var total = 0.02586            for value in data {2587                let u = (point - value) / bandwidth2588                if abs(u) < 1 {2589                    total += 0.75 * (1 - u * u)   // Epanechnikov kernel2590                }2591            }2592            densities[g] = total / (Double(n) * bandwidth)2593        }2594        return ZQPlotSpec.Series(label: label, x: xs, y: densities)2595    }25962597    // MARK: - Qualifiers25982599    /// Combined if/in keep-mask over current observations.2600    private func observationMask(_ command: ZQCommand) throws -> [Bool] {2601        var mask = [Bool](repeating: true, count: frame.rowCount)2602        if let condition = command.condition {2603            let evaluator = ExpressionEvaluator(frame: frame)2604            mask = try evaluator.evaluateCondition(condition)2605        }2606        if let range = command.range {2607            let n = frame.rowCount2608            func resolve(_ bound: Int) -> Int {2609                bound >= 0 ? bound : n + bound + 12610            }2611            let lower = max(1, resolve(range.lower))2612            let upper = min(n, resolve(range.upper))2613            for i in 0..<n where i + 1 < lower || i + 1 > upper {2614                mask[i] = false2615            }2616        }2617        return mask2618    }2619}2620