I know my answer is late, but might help others.
To answer your first question, yes each time series could be studied independently as univariate time series, by obtaining the mean and the autocovariance function for each series. However this approach doesn't take into account the possible dependence between the series.
By doing a linear regression to each series you are estimating the trend of the series, which does come into play when trying to fit a model for forecasting.
I will outline a very general approach for univariate time series that can be extended to multivariate time series:
Plot the series: check for trend and seasonality, changes in behaviour, outliers, etc.
Estimate the trend: a) with a smoothing procedure such as moving averages (no estimates) or b) model the trend with a regression equation.
"De-trend" the series. For additive models subtract the trend. For multiplicative models divide the series by the trend values.
Determine seasonal factors. The usual method is to average the "de-trended" values for a specific season.
Determine the random (residuals) component:
For an additive model: random = series - trend
For a multiplicative model: random = series/(trend*seasonal).
Choose a model to fit the residuals, using sample autocorrelation function.
Use residuals to forecast and then invert the transformations described above to arrive at forecasts of the original series.
You can check out Rob J Hyndman Forecasting Principles and Practice and the
The Little Book of Time Series for a better exposition.
With respect to your second question, you might want to use multivariate time series, in which a vector / matrix approach is used to mostly fit vector autoregresive models (VAR). You can find a much better explanation in Vector Autoregresive Models for Multivariate Time series, and how to use the R package vars in VAR, SVAR, SVEC Models: Implementation within R Package vars