Ruben the best way to proceed is to create a vector-valued ARIMA process with two components. Then based on this, you can build graphs as follows (replace "rFmode1" with the name of the process you defime on your own:

nssvmean = TimeSeriesThread[Mean, rFmode1["PathComponent", 1]];
spsvmean = TimeSeriesThread[Mean, rFmode1["PathComponent", 2]];
nsmup = TimeSeriesThread[Quantile[#, 0.95] &, rFmode1["PathComponent", 1]];
nsmdn = TimeSeriesThread[Quantile[#, 0.05] &, rFmode1["PathComponent", 1]];
spmup = TimeSeriesThread[Quantile[#, 0.95] &, rFmode1["PathComponent", 2]];
spmdn = TimeSeriesThread[Quantile[#, 0.05] &, rFmode1["PathComponent", 2]];

{DateListPlot[{indices["PathComponent", 1], 
   TimeSeries[nssvmean["Values"], {{2014, 11, 28}, Automatic}], 
   TimeSeries[nsmup["Values"], {{2014, 11, 28}, Automatic}], 
   TimeSeries[nsmdn["Values"], {{2014, 11, 28}, Automatic}]}, 
  Filling -> {3 -> {4}},
  PlotLabel -> "Nasdaq series with 3M projections bands"], 
 DateListPlot[{indices["PathComponent", 2], 
   TimeSeries[spsvmean["Values"], {{2014, 11, 28}, Automatic}], 
   TimeSeries[spmup["Values"], {{2014, 11, 28}, Automatic}], 
   TimeSeries[spmdn["Values"], {{2014, 11, 28}, Automatic}]}, 
  Filling -> {3 -> {4}},
  PlotLabel -> "S&P series with 3M projections bands"]}
DateListPlot[{indices["PathComponent", 2], 
  TimeSeries[spsvmean["Values"], {{2014, 11, 28}, Automatic}], 
  TimeSeries[spmup["Values"], {{2014, 11, 28}, Automatic}], 
  TimeSeries[spmdn["Values"], {{2014, 11, 28}, Automatic}]}, 
 Filling -> {3 -> {4}},
 PlotLabel -> "S&P series with 3M projections bands"]

{ListPlot[TimeSeriesThread[First[#] - Last[#] &, {nssvmean, nsmean}], 
  PlotStyle -> Blue, PlotLabel -> "Stoch CoVar Nasdaq mean difference"], 
 ListPlot[TimeSeriesThread[First[#] - Last[#] &, {spsvmean, spmean}], 
  PlotStyle -> Magenta, PlotLabel -> "Sotch CoVar S&P mean difference"]}

spath = RandomFunction[estPmod2, {1, 90}];
sts = TimeSeries[spath, ResamplingMethod -> Automatic];

sfunc = sts["PathFunction"];
gfunc[t_] := Join[{t}, sfunc[t]];

ParametricPlot3D[gfunc[t], {t, 1, 90}, 
 ColorFunction -> 
  Function[{t, x, y}, (ColorData["TemperatureMap"][t])], 
 PlotStyle -> Thickness[0.012], PlotTheme -> "Marketing", 
 BoxRatios -> 1, AxesLabel -> {t, "Nasdaq", "S&P"}]

Hi Igor, I'm going through your file trying to replicate the results in Mathematica but I'm having difficulties recreating all the graphs except for the first 3. Would it be possible to know the Mathematica code to generate the rest of the graphs? Many thanks in advance, Ruben