Retorno de un Activo Financiero

Según Fabozzi et al. (1996), un activo es cualquier posesión que tiene valor en un intercambio. Los activos pueden clasificarse como tangibles o intangibles.

Un activo tangible es aquel cuyo valor depende de sus propias características (edificios, maquinaria, terrenos, etc.).
Un activo intangible representa obligaciones legales sobre algún beneficio futuro.

Los activos financieros son activos intangibles, dado que su valor corresponde a una obligación monetaria futura.

Por ejemplo, un bono es un certificado de deuda, es decir, una promesa de pago futura que especifica monto, plazo, moneda y secuencia de pagos.
Cuando un inversionista compra un bono, está prestando dinero al emisor, quien se compromete a pagar intereses periódicos y devolver el capital al vencimiento.

Otro ejemplo son las acciones, que otorgan derecho a dividendos.
A diferencia de los bonos, el poseedor de una acción común no tiene certeza sobre el monto ni el momento en que se pagarán los dividendos.

Sea \(P_t\), para \(t = 1,2,\dots\), el precio de un activo en el tiempo \(t\).
Si el activo no paga dividendos, el retorno simple desde \(t-1\) hasta \(t\) se define como:

\[ R_t = \frac{P_t}{P_{t-1}} - 1 \]

o equivalentemente,

\[ R_t = \frac{P_t - P_{t-1}}{P_{t-1}} \]

El log-retorno (\(r_t\)) se define como:

\[ r_t = \ln(1 + R_t) \]

lo que equivale a:

\[ r_t = \ln(P_t) - \ln(P_{t-1}) \]

La serie de log-retornos no tiene unidades, es aproximadamente estable en media y facilita el cálculo de retornos compuestos en \(k\) períodos:

\[ r_t^{[k]} = \ln(1 + R_t^{[k]}) \]

y puede expresarse como la suma de log-retornos individuales:

\[ r_t^{[k]} = \sum_{i=0}^{k-1} r_{t-i} \]

Características de la Volatilidad

La volatilidad se define como la varianza condicional de la serie de retornos:

\[ \sigma_t^2 = \text{Var}(r_t \mid \mathcal{F}_{t-1}) \]

donde \(\mathcal{F}_{t-1}\) representa la información disponible hasta el período \(t-1\).

En series financieras, aunque los retornos puedan ser estacionarios con varianza incondicional constante, presentan variaciones de corto plazo en su varianza condicional, fenómeno conocido como heterocedasticidad.

Un buen modelo de volatilidad debe capturar:

Estos modelos son fundamentales para la gestión de riesgo, selección de portafolio y valoración de activos financieros.

# ===============================
# EJEMPLO: RETORNOS Y VOLATILIDAD
# ===============================

# Instalar paquetes si es necesario
# install.packages(c("quantmod","rugarch"))

library(quantmod)
library(rugarch)

# 1. Descargar datos (Apple)
getSymbols("AAPL", src = "yahoo", from = "2020-01-01")
## [1] "AAPL"
# Usamos precios de cierre ajustado
precio <- Ad(AAPL)

# 2. Calcular retornos simples
ret_simple <- dailyReturn(precio, type = "arithmetic")

# 3. Calcular log-retornos
ret_log <- dailyReturn(precio, type = "log")

# 4. Graficar precios y retornos
par(mfrow = c(3,1))

plot(precio, main="Precio Ajustado AAPL")
plot(ret_simple, main="Retorno Simple")
plot(ret_log, main="Log-Retorno")

par(mfrow = c(1,1))

# 5. Ajustar modelo GARCH(1,1)
spec <- ugarchspec(
  variance.model = list(model="sGARCH", garchOrder=c(1,1)),
  mean.model = list(armaOrder=c(0,0)),
  distribution.model = "norm"
)

modelo_garch <- ugarchfit(spec = spec, data = ret_log)

# 6. Mostrar resultados
show(modelo_garch)
## 
## *---------------------------------*
## *          GARCH Model Fit        *
## *---------------------------------*
## 
## Conditional Variance Dynamics    
## -----------------------------------
## GARCH Model  : sGARCH(1,1)
## Mean Model   : ARFIMA(0,0,0)
## Distribution : norm 
## 
## Optimal Parameters
## ------------------------------------
##         Estimate  Std. Error  t value Pr(>|t|)
## mu      0.001227    0.000411   2.9817 0.002866
## omega   0.000014    0.000001   9.9591 0.000000
## alpha1  0.096002    0.007515  12.7753 0.000000
## beta1   0.868144    0.011073  78.4005 0.000000
## 
## Robust Standard Errors:
##         Estimate  Std. Error  t value Pr(>|t|)
## mu      0.001227    0.000450   2.7253 0.006423
## omega   0.000014    0.000003   4.9344 0.000001
## alpha1  0.096002    0.011021   8.7107 0.000000
## beta1   0.868144    0.014677  59.1484 0.000000
## 
## LogLikelihood : 4044.382 
## 
## Information Criteria
## ------------------------------------
##                     
## Akaike       -5.2201
## Bayes        -5.2063
## Shibata      -5.2201
## Hannan-Quinn -5.2150
## 
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                     0.3371  0.5615
## Lag[2*(p+q)+(p+q)-1][2]    0.3385  0.7747
## Lag[4*(p+q)+(p+q)-1][5]    0.8506  0.8925
## d.o.f=0
## H0 : No serial correlation
## 
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                     0.3676  0.5443
## Lag[2*(p+q)+(p+q)-1][5]    2.2916  0.5513
## Lag[4*(p+q)+(p+q)-1][9]    4.3195  0.5356
## d.o.f=2
## 
## Weighted ARCH LM Tests
## ------------------------------------
##             Statistic Shape Scale P-Value
## ARCH Lag[3]    0.3006 0.500 2.000  0.5835
## ARCH Lag[5]    3.7738 1.440 1.667  0.1955
## ARCH Lag[7]    5.0502 2.315 1.543  0.2195
## 
## Nyblom stability test
## ------------------------------------
## Joint Statistic:  10.2755
## Individual Statistics:             
## mu     0.1634
## omega  1.9221
## alpha1 0.3216
## beta1  0.1598
## 
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic:          1.07 1.24 1.6
## Individual Statistic:     0.35 0.47 0.75
## 
## Sign Bias Test
## ------------------------------------
##                    t-value   prob sig
## Sign Bias           0.1624 0.8710    
## Negative Sign Bias  1.2415 0.2146    
## Positive Sign Bias  0.1225 0.9025    
## Joint Effect        2.0461 0.5629    
## 
## 
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
##   group statistic p-value(g-1)
## 1    20     69.16    1.269e-07
## 2    30     83.28    3.798e-07
## 3    40     83.78    4.107e-05
## 4    50    105.23    5.516e-06
## 
## 
## Elapsed time : 0.196574
# 7. Extraer volatilidad estimada
volatilidad <- sigma(modelo_garch)

# Graficar volatilidad condicional
plot(volatilidad, main="Volatilidad Condicional GARCH(1,1)")

Información de Fisher y Método de Newton–Raphson

Sea una muestra i.i.d. \(X_1,\dots,X_n\) con función de densidad \(f(x;\theta)\), donde \(\theta \in \mathbb{R}\) es un parámetro desconocido.

La función de verosimilitud es:

\[ L(\theta) = \prod_{i=1}^{n} f(X_i;\theta) \]

y el logaritmo de la verosimilitud:

\[ \ell(\theta) = \sum_{i=1}^{n} \log f(X_i;\theta) \]

1. Función Score

La función score se define como la derivada del log-verosimilitud:

\[ U(\theta) = \frac{\partial \ell(\theta)}{\partial \theta} \]

Bajo condiciones regulares,

\[ \mathbb{E}[U(\theta)] = 0 \]

2. Información de Fisher

La información de Fisher se define como:

\[ \mathcal{I}(\theta) = \mathbb{E}\left[ \left( \frac{\partial \ell(\theta)}{\partial \theta} \right)^2 \right] \]

Equivalentemente,

\[ \mathcal{I}(\theta) = - \mathbb{E}\left[ \frac{\partial^2 \ell(\theta)}{\partial \theta^2} \right] \]

La segunda igualdad se obtiene usando que:

\[ \frac{\partial}{\partial \theta} \mathbb{E}[U(\theta)] = 0 \]

y bajo intercambio de derivada e integral.

Para una muestra i.i.d.:

\[ \mathcal{I}_n(\theta) = n \mathcal{I}_1(\theta) \]

3. Expansión de Taylor

Para encontrar el estimador de máxima verosimilitud (MLE), resolvemos:

\[ U(\theta) = 0 \]

Aplicamos expansión de Taylor alrededor de \(\theta_0\):

\[ U(\theta) \approx U(\theta_0) + (\theta - \theta_0) U'(\theta_0) \]

Igualando a cero:

\[ 0 \approx U(\theta_0) + (\theta - \theta_0) U'(\theta_0) \]

Despejando:

\[ \theta \approx \theta_0 - \frac{U(\theta_0)}{U'(\theta_0)} \]

Como

\[ U'(\theta) = \frac{\partial^2 \ell(\theta)}{\partial \theta^2} \]

obtenemos el método de Newton–Raphson:

\[ \theta^{(k+1)} = \theta^{(k)} - \frac{U(\theta^{(k)})} {\frac{\partial^2 \ell(\theta^{(k)})}{\partial \theta^2}} \]

4. Método de Fisher Scoring

Usando que:

\[ \mathbb{E}\left[ \frac{\partial^2 \ell(\theta)}{\partial \theta^2} \right] = - \mathcal{I}(\theta) \]

reemplazamos el Hessiano por su esperanza:

\[ \theta^{(k+1)} = \theta^{(k)} + \frac{U(\theta^{(k)})} {\mathcal{I}(\theta^{(k)})} \]

Este procedimiento se conoce como Fisher Scoring.

5. Caso Multivariado

Si \(\theta \in \mathbb{R}^p\):

Score:

\[ U(\theta) = \nabla \ell(\theta) \]

Información de Fisher:

\[ \mathcal{I}(\theta) = - \mathbb{E} \left[ \nabla^2 \ell(\theta) \right] \]

Newton–Raphson:

\[ \theta^{(k+1)} = \theta^{(k)} - H^{-1}(\theta^{(k)}) U(\theta^{(k)}) \]

Fisher Scoring:

\[ \theta^{(k+1)} = \theta^{(k)} + \mathcal{I}^{-1}(\theta^{(k)}) U(\theta^{(k)}) \]

donde \(H(\theta)\) es el Hessiano observado.

6. Interpretación

  • Newton–Raphson usa el Hessiano observado.
  • Fisher Scoring usa la información esperada.
  • Ambos son métodos iterativos para resolver \(U(\theta)=0\).
  • Bajo regularidad, el MLE satisface:

\[ \sqrt{n}(\hat{\theta}-\theta_0) \overset{d}{\to} \mathcal{N} \left( 0,\mathcal{I}^{-1}(\theta_0) \right) \]

lo que muestra la conexión entre información de Fisher, eficiencia y varianza asintótica.

Modelos de Volatilidad Condicional

Sea \(r_t\) la serie de retornos financieros definida como:

\[ r_t = \mu + \varepsilon_t \]

donde:

\[ \varepsilon_t = \sigma_t z_t \]

con:

El objetivo es modelar la varianza condicional \(\sigma_t^2\).


1. Modelo ARCH(q)

Propuesto por Engle (1982).

Definición matemática

\[ \sigma_t^2 = \omega + \sum_{i=1}^{q} \alpha_i \varepsilon_{t-i}^2 \]

con restricciones:

\[ \omega > 0, \quad \alpha_i \ge 0 \]

Interpretación

  • La volatilidad depende solo de choques pasados.
  • Grandes errores pasados generan alta volatilidad futura.

Beneficios

  • Captura heterocedasticidad condicional.
  • Modelo simple y fácil de estimar.
  • Útil cuando la persistencia es baja.

2. Modelo GARCH(p,q)

Propuesto por Bollerslev (1986).

Definición matemática

\[ \sigma_t^2 = \omega + \sum_{i=1}^{q} \alpha_i \varepsilon_{t-i}^2 + \sum_{j=1}^{p} \beta_j \sigma_{t-j}^2 \]

Caso más usado: GARCH(1,1)

\[ \sigma_t^2 = \omega + \alpha \varepsilon_{t-1}^2 + \beta \sigma_{t-1}^2 \]

Restricciones:

\[ \omega > 0, \quad \alpha \ge 0, \quad \beta \ge 0 \]

Condición de estacionariedad:

\[ \alpha + \beta < 1 \]

Interpretación

  • Combina choques pasados y volatilidad pasada.
  • Captura persistencia en la volatilidad.

Beneficios

  • Reproduce aglomeración de volatilidad.
  • Modela memoria larga de forma parsimoniosa.
  • Es el modelo estándar en finanzas.

3. Modelo EGARCH(p,q)

Propuesto por Nelson (1991).

Definición matemática

Modela el logaritmo de la varianza:

\[ \ln(\sigma_t^2) = \omega + \sum_{i=1}^{q} \left[ \alpha_i \frac{\varepsilon_{t-i}}{\sigma_{t-i}} + \gamma_i \left( \left| \frac{\varepsilon_{t-i}}{\sigma_{t-i}} \right| - \mathbb{E}|z_t| \right) \right] + \sum_{j=1}^{p} \beta_j \ln(\sigma_{t-j}^2) \]

Característica clave

No requiere restricciones de positividad, ya que:

\[ \sigma_t^2 = \exp(\cdot) > 0 \]

Interpretación

  • Permite efectos asimétricos.
  • Si \(\alpha_i < 0\), malas noticias generan mayor volatilidad que buenas noticias (efecto leverage).

Beneficios

  • Captura asimetría.
  • Garantiza positividad automáticamente.
  • Más realista para mercados financieros.

Comparación Teórica

Modelo Persistencia Asimetría Restricciones
ARCH Baja No
GARCH Alta No
EGARCH Alta No

Resumen Conceptual

Estos modelos son fundamentales para:

# ==========================================
# EJEMPLO: ARCH vs GARCH vs EGARCH
# Evaluación de desempeño y bondad de ajuste
# ==========================================

# install.packages(c("quantmod","rugarch","tseries","FinTS"))
library(quantmod)
library(rugarch)
library(tseries)
library(FinTS)

# 1. Descargar datos (S&P 500)
getSymbols("^GSPC", src="yahoo", from="2018-01-01")
## [1] "GSPC"
precio <- Ad(GSPC)
ret_log <- dailyReturn(precio, type="log")
ret_log <- na.omit(ret_log)

# Separar muestra entrenamiento y prueba
n <- length(ret_log)
train <- ret_log[1:round(0.8*n)]
test  <- ret_log[(round(0.8*n)+1):n]

# =====================
# 2. ESPECIFICACIONES
# =====================

# ARCH(1)
spec_arch <- ugarchspec(
  variance.model=list(model="sGARCH", garchOrder=c(1,0)),
  mean.model=list(armaOrder=c(0,0)),
  distribution.model="norm"
)

# GARCH(1,1)
spec_garch <- ugarchspec(
  variance.model=list(model="sGARCH", garchOrder=c(1,1)),
  mean.model=list(armaOrder=c(0,0)),
  distribution.model="norm"
)

# EGARCH(1,1)
spec_egarch <- ugarchspec(
  variance.model=list(model="eGARCH", garchOrder=c(1,1)),
  mean.model=list(armaOrder=c(0,0)),
  distribution.model="norm"
)

# =====================
# 3. ESTIMACIÓN
# =====================

fit_arch  <- ugarchfit(spec_arch,  data=train)
fit_garch <- ugarchfit(spec_garch, data=train)
fit_egarch<- ugarchfit(spec_egarch,data=train)

# =====================
# 4. CRITERIOS DE INFORMACIÓN
# =====================

info <- rbind(
  ARCH  = infocriteria(fit_arch),
  GARCH = infocriteria(fit_garch),
  EGARCH= infocriteria(fit_egarch)
)

print(info)
##                       
## Akaike       -2.787315
## Bayes        -2.777437
## Shibata      -2.787322
## Hannan-Quinn -2.783652
## Akaike       -6.416619
## Bayes        -6.403449
## Shibata      -6.416631
## Hannan-Quinn -6.411734
## Akaike       -6.440850
## Bayes        -6.424387
## Shibata      -6.440868
## Hannan-Quinn -6.434745
# =====================
# 5. DIAGNÓSTICOS
# =====================

res_garch <- residuals(fit_garch, standardize=TRUE)

# Prueba Ljung-Box para autocorrelación
Box.test(res_garch, lag=10, type="Ljung-Box")
## 
##  Box-Ljung test
## 
## data:  res_garch
## X-squared = 10.599, df = 10, p-value = 0.3896
# Prueba ARCH LM
ArchTest(res_garch, lags=10)
## 
##  ARCH LM-test; Null hypothesis: no ARCH effects
## 
## data:  res_garch
## Chi-squared = 5.7174, df = 10, p-value = 0.8384
# =====================
# 6. PRONÓSTICO
# =====================

forecast_garch <- ugarchforecast(fit_garch, n.ahead=length(test))
vol_forecast <- sigma(forecast_garch)

# Error cuadrático medio en volatilidad (proxy: retornos^2)
realized_vol <- test^2
mse <- mean((realized_vol - vol_forecast^2)^2)

print(paste("MSE:", mse))
## [1] "MSE: 2.35234756507107e-07"
# =====================
# 7. GRÁFICAS
# =====================

par(mfrow=c(2,1))
plot(sigma(fit_garch), main="Volatilidad Condicional GARCH(1,1)")
plot(res_garch, main="Residuos Estandarizados GARCH")

par(mfrow=c(1,1))
# =========================
# 3. ESTIMACIÓN
# =========================

fit_arch   <- ugarchfit(spec_arch, data=train)
fit_garch  <- ugarchfit(spec_garch, data=train)
fit_egarch <- ugarchfit(spec_egarch, data=train)

# =========================
# 4. COMPARAR AIC y BIC
# =========================

criterios <- rbind(
  ARCH   = infocriteria(fit_arch),
  GARCH  = infocriteria(fit_garch),
  EGARCH = infocriteria(fit_egarch)
)

print(round(criterios,4))
##                     
## Akaike       -2.7873
## Bayes        -2.7774
## Shibata      -2.7873
## Hannan-Quinn -2.7837
## Akaike       -6.4166
## Bayes        -6.4034
## Shibata      -6.4166
## Hannan-Quinn -6.4117
## Akaike       -6.4409
## Bayes        -6.4244
## Shibata      -6.4409
## Hannan-Quinn -6.4347
# =========================
# 5. DIAGNÓSTICO RESIDUOS
# =========================

res_std <- residuals(fit_garch, standardize=TRUE)

Box.test(res_std, lag=20, type="Ljung-Box")
## 
##  Box-Ljung test
## 
## data:  res_std
## X-squared = 19.001, df = 20, p-value = 0.5217
Box.test(res_std^2, lag=20, type="Ljung-Box")
## 
##  Box-Ljung test
## 
## data:  res_std^2
## X-squared = 10.777, df = 20, p-value = 0.9518
# =========================
# 6. PRONÓSTICO
# =========================

forecast <- ugarchforecast(fit_garch, n.ahead=length(test))
vol_pred <- sigma(forecast)

# Proxy volatilidad realizada
realized_vol <- test^2

# MSE
mse <- mean((realized_vol - vol_pred^2)^2)
print(paste("MSE pronóstico:", round(mse,8)))
## [1] "MSE pronóstico: 2.4e-07"
# =========================
# 7. GRÁFICAS
# =========================

par(mfrow=c(2,1))
plot(vol_pred, type="l", col="blue",
     main="Volatilidad Pronosticada GARCH(1,1)")
plot(realized_vol, type="l", col="red",
     main="Volatilidad Realizada (r_t^2)")

par(mfrow=c(1,1))
# ==========================================================
# ESTUDIO COMPLETO: ARCH vs GARCH vs EGARCH
# Simulación + Estimación + Diagnóstico + Pronóstico
# ==========================================================

library(rugarch)
set.seed(123)

# ----------------------------------------------------------
# 1. SIMULACIÓN DEL PROCESO VERDADERO (GARCH(1,1))
# ----------------------------------------------------------

# Parámetros verdaderos
omega  <- 0.00001
alpha1 <- 0.08
beta1  <- 0.90
mu     <- 0

# Verificar estacionariedad
alpha1 + beta1   # Debe ser < 1
## [1] 0.98
spec_sim <- ugarchspec(
  variance.model = list(model="sGARCH", garchOrder=c(1,1)),
  mean.model     = list(armaOrder=c(0,0), include.mean=TRUE),
  distribution.model = "norm",
  fixed.pars = list(
    mu=mu,
    omega=omega,
    alpha1=alpha1,
    beta1=beta1
  )
)

sim <- ugarchpath(spec_sim, n.sim=2500)

ret <- fitted(sim)
sigma_true <- sigma(sim)

# Graficar serie simulada
par(mfrow=c(2,1))
plot(ret, type="l", main="Retornos simulados")
plot(sigma_true, type="l", main="Volatilidad verdadera")

par(mfrow=c(1,1))


# ----------------------------------------------------------
# 2. DIVISIÓN TRAIN / TEST
# ----------------------------------------------------------

train <- ret[1:2000]
test  <- ret[2001:2500]


# ----------------------------------------------------------
# 3. ESPECIFICACIÓN DE MODELOS
# ----------------------------------------------------------

# ARCH(1)
spec_arch <- ugarchspec(
  variance.model=list(model="sGARCH", garchOrder=c(1,0)),
  mean.model=list(armaOrder=c(0,0), include.mean=TRUE),
  distribution.model="norm"
)

# GARCH(1,1)
spec_garch <- ugarchspec(
  variance.model=list(model="sGARCH", garchOrder=c(1,1)),
  mean.model=list(armaOrder=c(0,0), include.mean=TRUE),
  distribution.model="norm"
)

# EGARCH(1,1)
spec_egarch <- ugarchspec(
  variance.model=list(model="eGARCH", garchOrder=c(1,1)),
  mean.model=list(armaOrder=c(0,0), include.mean=TRUE),
  distribution.model="norm"
)


# ----------------------------------------------------------
# 4. ESTIMACIÓN
# ----------------------------------------------------------

fit_arch   <- ugarchfit(spec_arch, data=train)
fit_garch  <- ugarchfit(spec_garch, data=train)
fit_egarch <- ugarchfit(spec_egarch, data=train)


# ----------------------------------------------------------
# 5. COMPARACIÓN DE CRITERIOS DE INFORMACIÓN
# ----------------------------------------------------------

info <- rbind(
  ARCH   = infocriteria(fit_arch),
  GARCH  = infocriteria(fit_garch),
  EGARCH = infocriteria(fit_egarch)
)

round(info,4)
##                     
## Akaike       -1.2220
## Bayes        -1.2136
## Shibata      -1.2220
## Hannan-Quinn -1.2189
## Akaike       -4.9059
## Bayes        -4.8947
## Shibata      -4.9059
## Hannan-Quinn -4.9018
## Akaike       -4.9010
## Bayes        -4.8870
## Shibata      -4.9010
## Hannan-Quinn -4.8959
# ----------------------------------------------------------
# 6. DIAGNÓSTICO DE RESIDUOS (GARCH)
# ----------------------------------------------------------

res_std <- residuals(fit_garch, standardize=TRUE)

# Autocorrelación
Box.test(res_std, lag=20, type="Ljung-Box")
## 
##  Box-Ljung test
## 
## data:  res_std
## X-squared = 17.31, df = 20, p-value = 0.6328
# Heterocedasticidad restante
Box.test(res_std^2, lag=20, type="Ljung-Box")
## 
##  Box-Ljung test
## 
## data:  res_std^2
## X-squared = 16.201, df = 20, p-value = 0.7041
# ----------------------------------------------------------
# 7. PRONÓSTICO OUT-OF-SAMPLE
# ----------------------------------------------------------

forecast_garch <- ugarchforecast(fit_garch, n.ahead=length(test))
vol_pred <- sigma(forecast_garch)

# Volatilidad realizada (proxy)
realized_vol <- test^2

# Error cuadrático medio
MSE <- mean((realized_vol - vol_pred^2)^2)
MSE
## [1] 2.402123e-07
# ----------------------------------------------------------
# 8. COMPARACIÓN VISUAL
# ----------------------------------------------------------

par(mfrow=c(2,1))
plot(vol_pred, type="l", col="blue",
     main="Volatilidad Pronosticada GARCH(1,1)")
plot(realized_vol, type="l", col="red",
     main="Volatilidad Realizada (r_t^2)")

par(mfrow=c(1,1))