Are the following steps correct to approach stationarity and cointegration?

Question

Say for 4 companies' close prices time-series with 500 observations each, I first estimate the Vector Autoregression of lagged order p, VAR(p). I achieved this by using Ordinary Linear Squares (OLS) with the following steps:

1) Form the dependent matrix Y with dimensions (row,columns): {4 x (500-p)}, where the oldest p observations are removed.

2) Form the explanatory matrix Z with dimensions: { (1+4*p) x (500-p )}

3) Now, calculate multivariate OLS estimator for regression coefficients matrix Beta as $\hat{B} = YZ'{(ZZ')}^{-1} $, where $Z'$ means transposed form of $Z$.

4) Then, estimate the regression residuals as: $\hat{res}=Y-\hat{B} *Z$

5) Determine the optimal lag via AIC and BIC (smallest value).

Now, I'm uncertain on how to go on from here. Do I repeat the same steps but this time by using the differenced time-series and then apply something like the Augmented-Dickey Fuller test on each individual time-series first i.e. test on company 1, 2, 3 and 4 separately? With this, determine if they are stationary or not by themselves. Is this the right way to apply the ADF test? I've read that ADF is applicable to a maximum of two time-series. Do I also test for the different combinations or would that be counter-productive?

Then, given that if they are non-stationary use the Johansen cointegration test to combine the differenced time-series to form a stationary process. (Even though this is unlikely, what happens if I get a mixture, i.e some are stationary or non-stationary, do I still use all of them?) Now, from what I understood about the Johansen framework, it requires a reduced rank matrix, usually defined as something like $m = \alpha\beta' $ which I believe is different from the echelon formed rank matrix I have in mind. For example in slide 23 of: http://statmath.wu.ac.at/~hauser/LVs/FinEtricsQF/FEtrics_Chp4.pdf

Am I right to say to obtain $m$, it is simply $m=-|I-\hat{B}|$, where $\hat{B}$ is the coefficient matrix at the optimal lagged period only, for example if optimal is $p=3$ I will only use $B_3$, $I$ is the identity matrix and $|...|$ means absolute value?

Update

The following was helpful in regards to a mixture of stationary and non-stationary time-series.

VAR or VECM for a mix of stationary and nonstationary variables

Answer 1

You should start from investigating the presence of unit roots in your time series. Assuming the maximum order of integration is 1, you may end up with three scenarios:

1. All have unit roots
2. Some have unit roots
3. None has a unit root

If 1., then investigate cointegration with the Johansen procedure. You may again end up with three scenarios:

(i) Cointegration of the full system
(ii) Cointegration of one or several subsystems
(iii) No cointegration

Follow "VAR or VECM for a mix of stationary and nonstationary variables" from there on.

If 2., also follow the post above.

If 3., estimate a VAR model for raw data.

Now, I'm uncertain on how to go on from here. Do I repeat the same steps but this time by using the differenced time-series and then apply something like the Augmented-Dickey Fuller test on each individual time-series first i.e. test on company 1, 2, 3 and 4 separately? With this, determine if they are stationary or not by themselves. Is this the right way to apply the ADF test? I've read that ADF is applicable to a maximum of two time-series. Do I also test for the different combinations or would that be counter-productive?

Regarding testing for unit roots in the individual series, you may use the ADF test. Regarding testing for cointegration, I recommend the Johansen procedure. Regarding ADF for more than one series, you are probably referring to the Engle-Granger procedure. If there is only one cointegrating vector, it will work (though you need to take the critical values from the appropriate null distribution which is not the standard one from the ADF test on a single series).

Then, given that if they are non-stationary use the Johansen cointegration test to combine the differenced time-series to form a stationary process. (Even though this is unlikely, what happens if I get a mixture, i.e some are stationary or non-stationary, do I still use all of them?)

See the post referenced above.

Now, from what I understood about the Johansen framework, it requires a reduced rank matrix, usually defined as something like $m = \alpha\beta' $ which I believe is different from the echelon formed rank matrix I have in mind. For example in slide 23 of: http://statmath.wu.ac.at/~hauser/LVs/FinEtricsQF/FEtrics_Chp4.pdf
Am I right to say to obtain $m$, it is simply $m=-|I-\hat{B}|$, where $\hat{B}$ is the coefficient matrix at the optimal lagged period only, for example if optimal is $p=3$ I will only use $B_3$, $I$ is the identity matrix and $|...|$ means absolute value?

I do not think $m=-|I-\hat{B}|$ is the relevant matrix. The relevant matrix comes out of the Johansen procedure. In general, you need to find the cointegrating vectors and then use their coefficients to weight the different cointegrated series to obtain stationary combinations thereof. E.g. in a bivariate system of cointegrated time series $(x,y)$, if the cointegrating vector is $(1,\beta)$, then a stationary combination of $x$ and $y$ will be $x+\beta y$.