title: "Shape Constraints with BsplineQuantReg"
author: "Alexandre Abbes"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Shape Constraints with BsplineQuantReg}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, fig.width = 7, fig.height = 5)
library(BsplineQuantReg)
```

## Introduction

This vignette covers shape constraints in B-spline quantile regression using the Karlin-Studden SOCP formulation. Shape constraints allow you to incorporate prior knowledge about the function's behavior (monotonicity, convexity, etc.) into the estimation process.

## Available Demos

The package includes several demos that illustrate shape constraints:

```{r demos}
# Example of available demos
dev.new()
demo(package = "BsplineQuantReg",temperature2)

# Run specific demos:
# demo("monotonicity")      - Increasing/decreasing constraints
# demo("convexity")         - Convexity/concavity constraints  
# demo("der3")              - Third derivative constraints
# demo("comprehensive")     - All constraint types
# demo("temperature")       - Real-world example with temperature data
```

## Uniform Constraints

### Monotonicity Constraints

Monotonicity constraints force the fitted function to be non-decreasing (`monot = 1`) or non-increasing (`monot = -1`).

```{r monotonicity}
# Generate data
set.seed(42)
n <- 100
x <- seq(0, 1, length.out = n)
y <- 2*x + 0.5*sin(4*pi*x) + 0.1*rnorm(n)
knots <- quantile(x, probs = seq(0, 1, length.out = 10))

# Unconstrained
fit_uncon <- quantile_spline(x, y, knots, tau = 0.5)

# Increasing constraint
fit_inc <- quantile_spline(x, y, knots, tau = 0.5, monot = 1)

# Decreasing constraint
fit_dec <- quantile_spline(x, y, knots, tau = 0.5, monot = -1)

# Compare
x_eval <- seq(0, 1, length.out = 200)
plot(x, y, pch = 16, cex = 0.5, col = "gray", main = "Monotonicity Constraints")
lines(x_eval, spline_eval(fit_uncon, x_eval), col = "red", lwd = 2, lty = 2)
lines(x_eval, spline_eval(fit_inc, x_eval), col = "blue", lwd = 2)
lines(x_eval, spline_eval(fit_dec, x_eval), col = "green", lwd = 2)
legend("topleft", legend = c("Unconstrained", "Increasing", "Decreasing"),
       col = c("red", "blue", "green"), lty = c(2, 1, 1), lwd = 2)
```

### Convexity Constraints

Convexity constraints force the second derivative to be non-negative (`convcons = 1`) or non-positive (`convcons = -1`).

```{r convexity}
# Convex function
y_conv <- (x - 0.5)^2 + 0.05*rnorm(n)

fit_convex <- quantile_spline(x, y_conv, knots, tau = 0.5, convcons = 1)

plot(x, y_conv, pch = 16, cex = 0.5, col = "gray", main = "Convexity Constraint")
lines(x_eval, spline_eval(fit_convex, x_eval), col = "blue", lwd = 2)
```

### Third Derivative Constraints

Third derivative constraints control the curvature of the second derivative. The control of its sign allows enforcing smoothness of the curvature.

```{r der3}
# Cubic function with varying third derivative
y_cubic <- x^3 - 1.5*x^2 + 0.5*x + 0.05*rnorm(n)

fit_der3_pos <- quantile_spline(x, y_cubic, knots, tau = 0.5, der3cons = 1)

plot(x, y_cubic, pch = 16, cex = 0.5, col = "gray", main = "Third Derivative Constraint")
lines(x_eval, spline_eval(fit_der3_pos, x_eval), col = "blue", lwd = 2)
```

### Constraints by Degree

| Degree | Monotonicity | Convexity | Third Derivative |
|--------|--------------|-----------|------------------|
| 1 (Linear) | +(linear) | X | X |
| 2 (Quadratic) | + (Karlin) | +(constant) | X|
| 3 (Cubic) | +(Karlin) | +(Karlin) | +(linear (constant)) |
| 4 (Quartic) | +(Karlin) | +(Karlin) | +(linear) |
'Karlin' means in fact 'quadratic' inequality.

## Partial Constraints

Partial constraints apply monotonicity or convexity only on specific intervals, not the entire domain.

### Partial Monotonicity

```{r partial-monot}
# Increasing only on first half
monot_partial <- c(rep(1, 5), rep(0, 4))  # 5 intervals increasing, 4 unconstrained

fit_partial <- quantile_spline(x, y, knots, tau = 0.5, monot = monot_partial)

plot(x, y, pch = 16, cex = 0.5, col = "gray", main = "Partial Monotonicity")
lines(x_eval, spline_eval(fit_partial, x_eval), col = "blue", lwd = 2)
abline(v = knots[6], col = "red", lty = 2)  # Boundary of constrained region
```

### Partial Convexity

```{r partial-convex}
# Convex on right half only
conv_partial <- rep(0, length(knots))
for (i in 1:length(knots)) {
  if (knots[i] > 0.5) conv_partial[i] <- 1
}

fit_conv_partial <- quantile_spline(x, y, knots, tau = 0.5, convcons = conv_partial)
plot(x, y, pch = 16, cex = 0.5, col = "gray", main = "Partial Convexity")
lines(x_eval, spline_eval(fit_conv_partial, x_eval), col = "blue", lwd = 2)
abline(v = knots[6], col = "red", lty = 2)  # Boundary of constrained region
```
## Mixed Constraints

Multiple shape constraints can be applied simultaneously.

```{r}
# Increasing AND convex
fit_both <- quantile_spline(x, y, knots, tau = 0.5, 
                            monot = 1, convcons = 1)

plot(x, y, pch = 16, cex = 0.5, col = "gray", 
     main = "Mixed Constraints: Increasing + Convex")
lines(x_eval, spline_eval(fit_both, x_eval), col = "blue", lwd = 2)
```

## Constraints with Multiple Knots

When using knot multiplicities, constraints need to be adapted to the reduced regularity. 

```{r multiple-knots}
# Create knots with multiplicity at a point
sn_mult <- c(0, 0, 0, 0, 0.3, 0.5, 0.5,0.5,0.5, 0.7, 1, 1, 1, 1)
basis_mult <- Bspline_base(sn_mult, degree = 3)
knots_mult <- basis_mult$knot

# Fit with constraints
fit_mult <- quantile_spline(x, y, knots_mult, tau = 0.5, monot = 1)
plot(x, y, pch = 16, cex = 0.5, col = "gray", 
     main = "Constraints with Multiple knots")
lines(x_eval, spline_eval(fit_mult, x_eval), col = "blue", lwd = 2)
```

## Visualizing Constraint Effects

```{r visualize-constraints}
par(mfrow = c(2, 2))

# No constraints
plot(x, y, pch = 16, cex = 0.4, col = "gray", main = "Unconstrained")
lines(x_eval, spline_eval(fit_uncon, x_eval), col = "red", lwd = 2)

# Monotonicity
plot(x, y, pch = 16, cex = 0.4, col = "gray", main = "Monotonicity")
lines(x_eval, spline_eval(fit_inc, x_eval), col = "blue", lwd = 2)

# Convexity
plot(x, y, pch = 16, cex = 0.4, col = "gray", main = "Convexity")
lines(x_eval, spline_eval(fit_convex, x_eval), col = "blue", lwd = 2)

# Mixed
plot(x, y, pch = 16, cex = 0.4, col = "gray", main = "Mixed")
lines(x_eval, spline_eval(fit_both, x_eval), col = "blue", lwd = 2)

par(mfrow = c(1, 1))
```
## Karlin-Studden SOCP Formulation

The shape constraints are implemented using the **Karlin-Studden characterization** of non-negative polynomials of degree 2 or 3. 

Depending on the degree of the spline (3 or 4), these constraints apply to the first or second derivative:

| Spline Degree | Constraint Type | Polynomial Degree | Applied to |
|---------------|-----------------|-------------------|------------|
| 3 (Cubic) | Monotonicity | Quadratic (degree 2) | First derivative |
| 3 (Cubic) | Convexity | Linear (degree 1) | Second derivative |
| 4 (Quartic) | Monotonicity | Cubic (degree 3) | First derivative |
| 4 (Quartic) | Convexity | Quadratic (degree 2) | Second derivative |

**Other constraints** (third derivative, constraints at knots) are **linear** and are added as linear inequalities.

The complete set of constraints forms a **Second-Order Cone Program (SOCP)**. This SOCP formulation ensures **exact shape constraints**, unlike approximations that only enforce constraints at knots or at a finite set of points.

The particular structure of the problem allows the use of the **CVXR** package, which handles this type of mixed quadratic/linear conic problems efficiently through its DCP (Disciplined Convex Programming) framework and supports multiple solvers (CLARABEL, OSQP, ECOS, SCS).

## Practical Considerations

### When to Use Shape Constraints

| Constraint | When to Use |
|------------|-------------|
| Monotonicity | Demand curves, growth curves, dose-response |
| Convexity | Cost functions, risk aversion, production functions |
| Third derivative | Smoothing curvature, spline regularity |
| Partial | Data with known local behavior |

### Constraint Strength

Constraints can be too strong if:
- The data clearly violates the assumed shape
- Too few knots leading to over-smoothing
- Multiple constraints conflicting

### Choosing the Degree

| Degree | Smoothness | Flexibility |
|--------|------------|-------------|
| 1 | Piecewise linear | Low |
| 2 | C¹ | Medium |
| 3 | C² | High |
| 4 | C³ | Very high |

## Summary

The package provides:
- **Uniform constraints**: Apply globally
- **Partial constraints**: Apply on specific intervals
- **Mixed constraints**: Multiple constraints simultaneously
- **SOCP formulation**: Exact (not approximated) constraints
- **multiple knots**: for lower regularity handling at knots (shocks)

For more examples, see the demos:

```{r}
demo(package = "BsplineQuantReg")
```
```

## Additional Sections You Could Add

| Section | Content |
|---------|---------|
| **Verification** | How to check if constraints are satisfied |
| **Comparison with cobs** | Differences in constraint implementation |
| **Real Examples** | Temperature, economics, biology applications |
| **Troubleshooting** | Common issues and solutions |