Introduction
It is pure in knowledge evaluation purposes for parameters to have bounds; variances cannot be damaging, GARCH coefficients should sum to lower than one for stationarity, and mixing proportions dwell between zero and one.
Whenever you estimate these fashions by most probability, the optimizer must respect these bounds, not simply on the answer, however all through the search. If optimization searches wander into invalid territory, it might impression the reliability and convergence of your outcomes. For instance, you might get complicated numbers from damaging variances, explosive forecasts from non-stationary GARCH, or likelihoods that make no sense.
GAUSS 26.0.1 introduces reduce, the primary new GAUSS optimizer in over 10 years, to deal with this cleanly.
The reduce optmizer let’s you specify bounds straight and GAUSS internally retains parameters possible at each iteration. No extra log-transforms, no penalty features, and no doublechecking.
In right now’s weblog, we’ll see the brand new reduce perform in motion, as we stroll via two examples:
- A GARCH estimation the place variance parameters should be optimistic
- A Stochastic frontier fashions the place each variance elements should be optimistic.
In each instances, bounded optimization makes estimation simpler and aligns outcomes with idea.
Why Bounds Matter
To see why this issues in follow, let’s have a look at a well-known instance. Think about a GARCH(1,1) mannequin:
$sigma^2_t = omega + alpha varepsilon^2_{t-1} + beta sigma^2_{t-1}$
For this mannequin to be well-defined and economically significant:
- The baseline variance should be optimistic ($omega gt 0$)
- Shocks and persistence should contribute non-negatively to variance ($alpha geq 0$, $beta geq 0$)
- The mannequin should be stationary ($alpha + beta lt 1$)
The normal workaround is to estimate remodeled parameters, $log(omega)$ as an alternative of $omega$, then convert again. This works, nevertheless it distorts the optimization floor and complicates customary error calculations. You are not estimating the parameters you care about; you are estimating transforms and hoping the numerics work out.
With bounded optimization, you estimate $omega$, $alpha$, and $beta$ straight, with the optimizer respecting the constraints all through.
Instance 1: GARCH(1,1) on Commodity Returns
Let’s estimate a GARCH(1,1) mannequin on a dataset of 248 observations of commodity value returns (this knowledge is included within the GAUSS 26 examples listing).
Step One: Information and Chance
First, we load the information and specify our log-likelihood goal perform.
// Load returns knowledge (ships with GAUSS)
fname = getGAUSShome("examples/df_returns.gdat");
returns = loadd(fname, "rcpi");
// GARCH(1,1) damaging log-likelihood
proc (1) = garch_negll(theta, y);
native omega, alpha, beta_, sigma2, ll, t;
omega = theta[1];
alpha = theta[2];
beta_ = theta[3];
sigma2 = zeros(rows(y), 1);
// Initialize with pattern variance
sigma2[1] = stdc(y)^2;
// Variance recursion
for t (2, rows(y), 1);
sigma2[t] = omega + alpha * y[t-1]^2 + beta_ * sigma2[t-1];
endfor;
// Gaussian log-likelihood
ll = -0.5 * sumc(ln(2*pi) + ln(sigma2) + (y.^2) ./ sigma2);
retp(-ll); // Return damaging for minimization
endp;
Step Two: Setting Up Optimization
Now we arrange the bounded optimization with:
- $omega gt 0$ (small optimistic decrease certain to keep away from numerical points)
- $alpha geq 0$
- $beta geq 0$
As a result of reduce handles easy field constraints, we impose particular person higher bounds on $alpha$ and $beta$ to maintain the optimizer in an inexpensive area. We’ll confirm the stationarity situation, $alpha + beta lt 1$ after estimation.
// Beginning values
theta0 = { 0.00001, // omega (small, let knowledge converse)
0.05, // alpha
0.90 }; // beta
// Arrange reduce
struct minimizeControl ctl;
ctl = minimizeControlCreate();
// Bounds: all parameters optimistic, alpha + beta < 1
ctl.bounds = { 1e-10 1, // omega in [1e-10, 1]
0 1, // alpha in [0, 1]
0 0.9999 }; // beta in [0, 0.9999]
We cap $beta$ barely under 1 to keep away from numerical points close to the boundary, the place the probability floor can turn into flat and unstable.
Step Three: Operating the Mannequin
Lastly, we name reduce to run our mannequin.
// Estimate
struct minimizeOut out;
out = reduce(&garch_negll, theta0, returns, ctl);
Outcomes and Visualization
After estimation, we’ll extract the conditional variance collection and ensure the stationarity situation:
// Extract estimates
omega_hat = out.x[1];
alpha_hat = out.x[2];
beta_hat = out.x[3];
print "omega = " omega_hat;
print "alpha = " alpha_hat;
print "beta = " beta_hat;
print "alpha + beta = " alpha_hat + beta_hat;
print "Iterations: " out.iterations;
Output:
omega = 0.0000070 alpha = 0.380 beta = 0.588 alpha + beta = 0.968 Iterations: 39
There are just a few noteworthy outcomes:
- The excessive persistence ($alpha + beta approx 0.97$) means volatility shocks decay slowly.
- The comparatively excessive $alpha$ (0.38) signifies that latest shocks have substantial fast impression on variance.
- The optimization converged in 39 iterations with all parameters staying inside their bounds all through. No invalid variance evaluations, no numerical exceptions.
Visualizing the conditional variance alongside the unique collection gives additional perception:
// Compute conditional variance collection for plotting
T = rows(returns);
sigma2_hat = zeros(T, 1);
sigma2_hat[1] = stdc(returns)^2;
for t (2, T, 1);
sigma2_hat[t] = omega_hat + alpha_hat * returns[t-1]^2 + beta_hat * sigma2_hat[t-1];
endfor;
// Plot returns and conditional volatility
struct plotControl plt;
plt = plotGetDefaults("xy");
plotSetTitle(&plt, "GARCH(1,1): Returns and Conditional Volatility");
plotSetYLabel(&plt, "Returns / Volatility");
plotLayout(2, 1, 1);
plotXY(plt, seqa(1, 1, T), returns);
plotLayout(2, 1, 2);
plotSetTitle(&plt, "Conditional Normal Deviation");
plotXY(plt, seqa(1, 1, T), sqrt(sigma2_hat));
The plot exhibits volatility clustering: durations of excessive volatility are likely to persist, in line with what we observe in commodity markets.
Instance 2: Stochastic Frontier Mannequin
Stochastic frontier evaluation separates random noise from systematic inefficiency. It is extensively utilized in productiveness evaluation to measure how far corporations function under their manufacturing frontier.
The mannequin:
$y = Xbeta + v – u$
the place:
- $v sim N(0, sigma^2_v)$ — symmetric noise (measurement error, luck)
- $u sim N^+(0, sigma^2_u)$ — one-sided inefficiency (all the time reduces output)
Each variance elements should be optimistic. If the optimizer tries $sigma^2_v lt 0$ or $sigma^2_u lt 0$, the probability includes sq. roots of damaging numbers.
Step One: Information and Chance
For this instance, we’ll simulate knowledge from a Cobb-Douglas manufacturing perform with inefficiency. This retains the instance self-contained and allows you to see precisely what’s being estimated.
// Simulate manufacturing knowledge
rndseed 8675309;
n = 500;
// Inputs (labor, capital, supplies)
labor = exp(2 + 0.5*rndn(n, 1));
capital = exp(3 + 0.7*rndn(n, 1));
supplies = exp(2.5 + 0.4*rndn(n, 1));
// True parameters
beta_true = { 1.5, // fixed
0.4, // labor elasticity
0.3, // capital elasticity
0.25 }; // supplies elasticity
sig2_v_true = 0.02; // noise variance
sig2_u_true = 0.08; // inefficiency variance
// Generate output with noise (v) and inefficiency (u)
v = sqrt(sig2_v_true) * rndn(n, 1);
u = sqrt(sig2_u_true) * abs(rndn(n, 1)); // half-normal
X = ones(n, 1) ~ ln(labor) ~ ln(capital) ~ ln(supplies);
y = X * beta_true + v - u; // inefficiency reduces output
After simulating our knowledge, we specify the log-likelihood perform for minimization:
// Stochastic frontier log-likelihood (half-normal inefficiency)
proc (1) = sf_negll(theta, y, X);
native ok, beta_, sig2_v, sig2_u, sigma, lambda;
native eps, z, ll;
ok = cols(X);
beta_ = theta[1:k];
sig2_v = theta[k+1];
sig2_u = theta[k+2];
sigma = sqrt(sig2_v + sig2_u);
lambda = sqrt(sig2_u / sig2_v);
eps = y - X * beta_;
z = -eps * lambda / sigma;
ll = -0.5*ln(2*pi) + ln(2) - ln(sigma)
- 0.5*(eps./sigma).^2 + ln(cdfn(z));
retp(-sumc(ll));
endp;
Step Two: Setting Up Optimization
As we did in our earlier instance, we start with our beginning values. For this mannequin, we run OLS and use the residual variance as beginning values:
// OLS for beginning values
beta_ols = invpd(X'X) * X'y;
resid = y - X * beta_ols;
sig2_ols = meanc(resid.^2);
// Beginning values: Cut up residual variance
// between noise and inefficiency
theta0 = beta_ols | (0.5 * sig2_ols) | (0.5 * sig2_ols);
We depart our coefficients unbounded however constrain the variances to be optimistic:
// Bounds: coefficients unbounded, variances optimistic
ok = cols(X);
struct minimizeControl ctl;
ctl = minimizeControlCreate();
ctl.bounds = (-1e300 * ones(ok, 1) | 0.001 | 0.001) ~ (1e300 * ones(ok+2, 1));
Step Three: Operating the Mannequin
Lastly, we name reduce to estimate our mannequin:
// Estimate
struct minimizeOut out;
out = reduce(&sf_negll, theta0, y, X, ctl);
Outcomes and Visualization
Now that we have estimated our mannequin, let’s look at our outcomes.
// Extract estimates
ok = cols(X);
beta_hat = out.x[1:k];
sig2_v_hat = out.x[k+1];
sig2_u_hat = out.x[k+2];
print "Coefficients:";
print " fixed = " beta_hat[1];
print " ln(labor) = " beta_hat[2];
print " ln(capital) = " beta_hat[3];
print " ln(supplies)= " beta_hat[4];
print "";
print "Variance elements:";
print " sig2_v (noise) = " sig2_v_hat;
print " sig2_u (inefficiency)= " sig2_u_hat;
print " ratio sig2_u/whole = " sig2_u_hat / (sig2_v_hat + sig2_u_hat);
print "";
print "Iterations: " out.iterations;
This prints out coefficients and variance elements:
Coefficients: fixed = 1.51 ln(labor) = 0.39 ln(capital) = 0.31 ln(supplies)= 0.24 Variance elements: sig2_v (noise) = 0.022 sig2_u (inefficiency)= 0.087 ratio sig2_u/whole = 0.80 Iterations: 38
The estimates recuperate the true parameters moderately nicely. The variance ratio ($approx 0.80$) tells us that the majority residual variation is systematic inefficiency, not measurement error — an vital discovering for coverage.
We are able to additionally compute and plot firm-level effectivity scores:
// Compute effectivity estimates (Jondrow et al. 1982)
eps = y - X * beta_hat;
sigma = sqrt(sig2_v_hat + sig2_u_hat);
lambda = sqrt(sig2_u_hat / sig2_v_hat);
mu_star = -eps * sig2_u_hat / (sig2_v_hat + sig2_u_hat);
sig_star = sqrt(sig2_v_hat * sig2_u_hat / (sig2_v_hat + sig2_u_hat));
// E[u|eps] - conditional imply of inefficiency
u_hat = mu_star + sig_star * (pdfn(mu_star/sig_star) ./ cdfn(mu_star/sig_star));
// Technical effectivity: TE = exp(-u)
TE = exp(-u_hat);
// Plot effectivity distribution
struct plotControl plt;
plt = plotGetDefaults("hist");
plotSetTitle(&plt, "Distribution of Technical Effectivity");
plotSetXLabel(&plt, "Technical Effectivity (1 = frontier)");
plotSetYLabel(&plt, "Frequency");
plotHist(plt, TE, 20);
print "Imply effectivity: " meanc(TE);
print "Min effectivity: " minc(TE);
print "Max effectivity: " maxc(TE);
Imply effectivity: 0.80 Min effectivity: 0.41 Max effectivity: 0.95
The histogram exhibits substantial variation in effectivity — some corporations function close to the frontier (TE $approx$ 0.95), whereas others produce 40-50% under their potential. That is the sort of perception that drives productiveness analysis.
Each variance estimates stayed optimistic all through optimization. No log-transforms wanted, and the estimates apply on to the parameters we care about.
When to Use reduce
The reduce process is designed for one factor: optimization with certain constraints. If that is all you want, it is the appropriate device.
| State of affairs | Suggestion |
|---|---|
| Parameters with easy bounds | reduce |
| Nonlinear constraints ($g(x) leq 0$) | sqpSolveMT |
| Equality constraints | sqpSolveMT |
| Algorithm switching, complicated issues | OPTMT |
For the GARCH and stochastic frontier examples above — and most MLE issues the place parameters have pure bounds — reduce handles it straight.
Conclusion
Bounded parameters present up continually in econometric fashions: variances, volatilities, possibilities, shares. GAUSS 26.0.1 provides you a clear solution to deal with them with reduce. As we noticed right now reduce:
- Set bounds within the management construction
- Optimizer respects bounds all through (not simply on the answer)
- No log-transforms or penalty features
- Included in base GAUSS
If you happen to’ve been working round parameter bounds with transforms or checking for invalid values inside your probability perform, that is the cleaner path.











