VAR - ARCH LM Test results are conflicting

I'm using the R package vars to model multivariate series with VAR of order p = 5. The multivariate series is: $$Y_t = \{y_{1,t}; y_{2,t}; y_{3,t}\}^\top$$

In order to stabilize the time series, I calculated the relative differences (return) of each element:

$$x_{i, t} = \frac{y_{t} - y_{t-1}}{y_{t-1}}$$

When I test the residuals for heteroskedasticity with the function arch.test, which implements an ARCH-LM test, I get conflicting results.

For example, with arch.test(data.var, lags.multi = 1), I get:

> ht <- arch.test(data.var, lags.multi = 1)
> ht

ARCH (multivariate)

data:  Residuals of VAR object data.var
Chi-squared = 27.663, df = 36, p-value = 0.839


which tells me to accept the null hypothesis. But when I run arch.test(data.var, lags.multi = 2), then I get:

> ht <- arch.test(data.var, lags.multi = 2)
> ht

ARCH (multivariate)

data:  Residuals of VAR object data.var
Chi-squared = 295.84, df = 72, p-value < 2.2e-16


now the result is that I should reject the null hypothesis and hence my multivariate time series is heteroskedastic. If I continue increasing the lags.multi parameter, the p-value closes to zero.

Why different test lags result in a different result? Should I use another method to stabilize my time series?

The ARCH-LM test (be it multivariate or univariate) with $q$ lags tests whether there are ARCH effects at lags from 1 up to $q$. It tests the joint significance of coefficients $\alpha_1,\dotsc,\alpha_q$ in the equation

$$x_t^2 = \alpha_0 + \alpha_1 x_{t-1}^2 + \dotsc + \alpha_q x_{t-q}^2 + \varepsilon_t$$

(I use different notation than in the Wikipedia article linked to above to stress that the dependent variable need not be a residual from an ARIMA model and to match your notation).

For $q=1$, it tests for ARCH effects at lag 1: $\text{H}_0:\alpha_1=0$.
For $q=2$, it tests for ARCH effects at lags 1 and 2 jointly: $\text{H}_0:\alpha_1=\alpha_2=0$.
These are different hypotheses so one could not say the result of the first test is in conflict with the result of the second test.

In your case, you do not find an ARCH effect at lag 1 but you do find them at lags 1 and 2 jointly. That suggests that the squared time series is autocorrelated at lag 2 but not lag 1.

If I continue increasing the lags.multi parameter, the p-value closes to zero.

That suggests the squared series has significant partial autocorrelations (significant ARCH effects) at higher lag orders as well; if that was not the case and the only significant autocorrelation was at lag 2, the $p$-values of the ARCH-LM test would decrease with increasing lag order as the single effect would be "diluted" in a test of joint significance of an increasing number of coefficients.

In sum, you seem to have an ARCH effect at lag 2 and likely higher lag orders, which means your data is conditionally heteroskedastic.