python - tsa - statsmodels time series

Intervalos de confianza y predicción con StatsModels.

Question

Hago esta regresión lineal con StatsModels:

import numpy as np
import statsmodels.api as sm
from statsmodels.sandbox.regression.predstd import wls_prediction_std

n = 100

x = np.linspace(0, 10, n)
e = np.random.normal(size=n)
y = 1 + 0.5*x + 2*e
X = sm.add_constant(x)

re = sm.OLS(y, X).fit()
print(re.summary())

prstd, iv_l, iv_u = wls_prediction_std(re)

Mis preguntas son: ¿ iv_l y iv_u son los intervalos de confianza superior o inferior o los intervalos de predicción? ¿Cómo consigo otros? Necesito la confianza y los intervalos de predicción de todos los puntos, para hacer una trama.

Answer 1

Para datos de prueba puede intentar usar lo siguiente.

predictions = result.get_prediction(out_of_sample_df)
predictions.summary_frame(alpha=0.05)

Encontré el método summary_frame () enterrado here y puede encontrar el método get_prediction () here . Puede cambiar el nivel de significación del intervalo de confianza y el intervalo de predicción modificando el parámetro "alfa".

Estoy publicando esto aquí porque esta fue la primera publicación que aparece cuando se busca una solución para intervalos de confianza y predicción, aunque esto concierne a los datos de prueba.

Aquí hay una función para tomar un modelo, nuevos datos y un cuantil arbitrario, utilizando este enfoque:

def ols_quantile(m, X, q):
  # m: OLS model.
  # X: X matrix.
  # q: Quantile.
  #
  # Set alpha based on q.
  a = q * 2
  if q > 0.5:
    a = 2 * (1 - q)
  predictions = m.get_prediction(X)
  frame = predictions.summary_frame(alpha=a)
  if q > 0.5:
    return frame.obs_ci_upper
  return frame.obs_ci_lower

Answer 2

Puede calcularlos basándose en los resultados dados por statsmodel y los supuestos de normalidad.

Aquí hay un ejemplo para OLS y CI para el valor medio:

import statsmodels.api as sm
import numpy as np
from scipy import stats

#Significance level:
sl = 0.05
#Evaluate mean value at a required point x0. Here, at the point (0.0,2.0) for N_model=2:
x0 = np.asarray([1.0, 0.0, 2.0])# If you have no constant in your model, remove the first 1.0. For more dimensions, add the desired values.

#Get an OLS model based on output y and the prepared vector X (as in your notation):
model = sm.OLS(endog = y, exog = X )
results = model.fit()
#Get two-tailed t-values:
(t_minus, t_plus) = stats.t.interval(alpha = (1.0 - sl), df =  len(results.resid) - len(x0) )
y_value_at_x0 = np.dot(results.params, x0)
lower_bound = y_value_at_x0 + t_minus*np.sqrt(results.mse_resid*( np.dot(np.dot(x0.T,results.normalized_cov_params),x0) ))
upper_bound = y_value_at_x0 +  t_plus*np.sqrt(results.mse_resid*( np.dot(np.dot(x0.T,results.normalized_cov_params),x0) ))

Puede envolver una buena función alrededor de esto con los resultados de entrada, el punto x0 y el nivel de significación sl.

Ahora no estoy seguro si puedes usar esto para WLS () ya que allí suceden cosas adicionales.

Ref: Ch3 en [DC Montgomery y EA Peck. “Introducción al Análisis de Regresión Lineal”. 4to. Ed., Wiley, 1992].

Answer 3

Puede obtener los intervalos de predicción utilizando la clase LRPI () del cuaderno Ipython en mi repositorio ( https://github.com/shahejokarian/regression-prediction-interval ).

Debe establecer el valor de t para obtener el intervalo de confianza deseado para los valores de predicción, de lo contrario, el valor predeterminado es 95% conf. intervalo.

La clase LRPI usa las bibliotecas LinearRegression, numpy y pandas de sklearn.linear_model.

Hay un ejemplo que se muestra en el cuaderno también.

Answer 4

iv_l y iv_u dan los límites del intervalo de predicción para cada punto.

El intervalo de predicción es el intervalo de confianza para una observación e incluye la estimación del error.

Creo que el intervalo de confianza para la predicción media aún no está disponible en statsmodels . (En realidad, el intervalo de confianza para los valores ajustados se esconde dentro de la tabla de sumario de influencia_utilidad, pero debo verificar esto).

Los métodos de predicción adecuados para los modelos de estadísticas están en la lista de tareas pendientes.

Adición

Los intervalos de confianza están disponibles para OLS, pero el acceso es un poco torpe.

Para ser incluido después de ejecutar su script:

from statsmodels.stats.outliers_influence import summary_table

st, data, ss2 = summary_table(re, alpha=0.05)

fittedvalues = data[:, 2]
predict_mean_se  = data[:, 3]
predict_mean_ci_low, predict_mean_ci_upp = data[:, 4:6].T
predict_ci_low, predict_ci_upp = data[:, 6:8].T

# Check we got the right things
print np.max(np.abs(re.fittedvalues - fittedvalues))
print np.max(np.abs(iv_l - predict_ci_low))
print np.max(np.abs(iv_u - predict_ci_upp))

plt.plot(x, y, ''o'')
plt.plot(x, fittedvalues, ''-'', lw=2)
plt.plot(x, predict_ci_low, ''r--'', lw=2)
plt.plot(x, predict_ci_upp, ''r--'', lw=2)
plt.plot(x, predict_mean_ci_low, ''r--'', lw=2)
plt.plot(x, predict_mean_ci_upp, ''r--'', lw=2)
plt.show()

Esto debería dar los mismos resultados que SAS, http://jpktd.blogspot.ca/2012/01/nice-thing-about-seeing-zeros.html

Answer 5

summary_frame y summary_table funcionan bien cuando necesita resultados exactos para un solo cuantil, pero no vectoriza bien. Esto proporcionará una aproximación normal del intervalo de predicción (no el intervalo de confianza) y funciona para un vector de cuantiles:

def ols_quantile(m, X, q):
  # m: Statsmodels OLS model.
  # X: X matrix of data to predict.
  # q: Quantile.
  #
  from scipy.stats import norm
  mean_pred = m.predict(X)
  se = np.sqrt(m.scale)
  return mean_pred + norm.ppf(q) * se