Checking heteroscedasticity in arima residuals

Question

I'm using auto.arima function to analyze my data. And here's what I get:

Series: data 
ARIMA(2,1,2) with drift         

Coefficients:
         ar1      ar2      ma1     ma2   drift
      1.6679  -0.8005  -1.2424  0.3125  0.4225
s.e.  0.1007   0.1107   0.1396  0.1413  0.1895

sigma^2 estimated as 17.34:  log likelihood=-438.37
AIC=888.73   AICc=889.3   BIC=906.99

I'm trying to perform a heteroscedastic test to my residual model using bptest from lmtest package with this code:

> bp<-bptest(lm(residuals(model)~1))
> bp

    studentized Breusch-Pagan test

data:  lm(residuals(model) ~ 1)
BP = 1.231e-30, df = 0, p-value < 2.2e-16

Am I doing bptest right? When I analyze another data with this code I always get the same df and p-value.

edit: here is the data :

    World_Oil_Prices
[1-10] [11-20] [21-30] [31-40] [41-50] [51-60] [61-70] [71-80] [81-90] [91-100] [101-110] [111-120] [121-130] [131-140] [141-150] [151-156]
17.79   22.25   18.73   12.72   16.12   27.49   25.95   18.69   28.28   28.59   37.63   48.75   59.67   53.53   75.91   131.22
17.69   23.51   20.12   12.49   16.24   23.45   27.24   18.52   27.53   29.68   35.54   46.00   54.17   57.22   81.27   121.87
19.46   23.29   19.16   13.80   18.75   27.23   25.02   19.15   24.79   26.88   37.93   43.67   56.63   50.14   90.54   96.85
20.78   20.54   17.24   13.26   20.21   29.62   25.66   19.98   27.89   29.01   42.08   52.55   66.85   54.46   89.76   69.16
19.12   19.42   15.07   11.88   22.37   28.16   27.55   23.64   30.77   29.12   41.65   52.24   63.49   57.78   85.53   46.03
18.56   17.98   14.18   10.41   22.19   29.41   26.97   25.43   32.88   29.95   46.87   58.74   62.26   63.25   93.51   38.60
19.56   19.47   13.24   11.32   24.22   32.08   24.80   25.69   30.36   31.40   42.23   58.20   68.08   66.75   99.32   
20.19   18.02   13.39   10.75   25.01   31.40   25.81   24.49   25.49   31.32   39.09   53.32   66.45   68.29   111.03  
22.14   18.45   13.97   12.86   25.21   32.33   25.03   25.75   26.06   33.67   42.76   49.41   56.38   73.69   123.35  
23.43   18.79   12.48   15.73   27.15   25.28   20.73   26.78   27.91   33.71   44.35   51.66   53.58   67.10   128.33  

And this is the R code:

library(forecast)
model<-auto.arima(data)
model
library(lmtest)
bp<-bptest(lm(residuals(model)~1))
bp

Answer 1

It appears to me that your model is possibly over-parameterized or over-differenced or both. This is the result of selecting from a set of fixed models and or the total ignoring of the effect of Pulses/Level Shifts/Seasonal Pulses/Local Time Trends. If you post your original data and your residual series I will test the constancy of the error variance by employing a search scheme along the lines suggested by TSAY http://www.unc.edu/~jbhill/tsay.pdf . Ignored by many researchers except a few since the mainstream time series software developers found it too difficult/technical to implement. Note that if the parameters of the model change over time this is often misdiagnosed as "evidence" of variance change. One needs to carefully diagnose possible Gauussian Violations.

Answer 2

After a quick look at the essentials of the BP test, it does not seem suitable to use a model of the form lm(resid~1) as an argument to the bptest function. The idea of the BP test is that error variance may vary together with an explanatory variable. That variable cannot be identically equal to a constant, otherwise the error variance could not vary with it. So having 1 on the right hand side of resid~1 does not make sense to me.

If you have an ARIMA model, you should perhaps follow the advice of @javlacalle here.
You should construct the variables corresponding to the AR(1), AR(2), MA(1) and MA(2) terms, say, $x_{t-1}$, $x_{t-2}$, $\varepsilon_{t-1}$, $\varepsilon_{t-2}$. Then build a model

$$\varepsilon_t^2=\beta_0+\beta_1 x_{t-1}+\beta_2 x_{t-2}+\beta_3 \varepsilon_{t-1}+\beta_4 \varepsilon_{t-2}+u_t$$

and test the model's overall significance (normally reported as $F$-statistic in model summary).
If the model is insignificant ($F$-statistic relatively low, associated $p$-value relatively high), there is not enough evidence to reject the null hypothesis "no conditional heteroskedasticity".
If the model is significant ($F$-statistic relatively high, associated $p$-value relatively low), the null hypothesis can be rejected in favour of the alternative hypothesis "presence of conditional heteroskedasticity".

I do not know the BP test very well and I do not remember ever seeing the BP test applied on residuals from an ARIMA model, so please correct me if I am wrong.