Makes forecast with the estimated model

        Parameters
        ----------
        h : int (default : 5)
            How many steps ahead would you like to forecast?

        past_values : int (default : 20)
            How many past observations to show on the forecast graph?

        intervals : Boolean
            Would you like to show prediction intervals for the forecast?

        oos_data : pd.DataFrame
            Data for the variables to be used out of sample (ys can be NaNs)

        Returns
        ----------
        - Plot of the forecast

def plot_predict(self, h=5, past_values=20, intervals=True, oos_data=None, **kwargs):
        """ Makes forecast with the estimated model

        Parameters
        ----------
        h : int (default : 5)
            How many steps ahead would you like to forecast?

        past_values : int (default : 20)
            How many past observations to show on the forecast graph?

        intervals : Boolean
            Would you like to show prediction intervals for the forecast?

        oos_data : pd.DataFrame
            Data for the variables to be used out of sample (ys can be NaNs)

        Returns
        ----------
        - Plot of the forecast
        """
        import matplotlib.pyplot as plt
        import seaborn as sns

        figsize = kwargs.get('figsize',(10,7))

        if self.latent_variables.estimated is False:
            raise Exception("No latent variables estimated!")
        else:
            # Sort/manipulate the out-of-sample data
            _, X_oos = dmatrices(self.formula, oos_data)
            X_oos = np.array([X_oos])[0]
            X_pred = X_oos[:h]

            date_index = self.shift_dates(h)

            if self.latent_variables.estimation_method in ['M-H']:
                
                sim_vector = np.zeros([15000,h])

                for n in range(0, 15000):
                    t_z = self.draw_latent_variables(nsims=1).T[0]
                    _, Y, _, coefficients = self._model(t_z)
                    coefficients_star = coefficients.T[-1]
                    theta_pred = np.dot(np.array([coefficients_star]), X_pred.T)[0]
                    t_z = np.array([self.latent_variables.z_list[k].prior.transform(t_z[k]) for k in range(t_z.shape[0])])
                    model_scale, model_shape, model_skewness = self._get_scale_and_shape(t_z)
                    sim_vector[n,:] = self.family.draw_variable(self.link(theta_pred), model_scale, model_shape, model_skewness, theta_pred.shape[0])
                mean_values = np.append(Y, self.link(np.array([np.mean(i) for i in sim_vector.T])))
            else:

                # Retrieve data, dates and (transformed) latent variables
                _, Y, _, coefficients = self._model(self.latent_variables.get_z_values()) 
                coefficients_star = coefficients.T[-1] 
                theta_pred = np.dot(np.array([coefficients_star]), X_pred.T)[0]  
                t_z = self.transform_z()
                sim_vector = np.zeros([15000,h])
                mean_values = np.append(Y, self.link(theta_pred))
                model_scale, model_shape, model_skewness = self._get_scale_and_shape(t_z)

                if self.model_name2 == "Skewt":
                    m1 = (np.sqrt(model_shape)*sp.gamma((model_shape-1.0)/2.0))/(np.sqrt(np.pi)*sp.gamma(model_shape/2.0))
                    mean_values += (model_skewness - (1.0/model_skewness))*model_scale*m1 

                for n in range(0,15000):
                    sim_vector[n,:] = self.family.draw_variable(self.link(theta_pred),model_scale,model_shape,model_skewness,theta_pred.shape[0])

            sim_vector = sim_vector.T
            error_bars = []
            for pre in range(5,100,5):
                error_bars.append(np.insert([np.percentile(i,pre) for i in sim_vector], 0, mean_values[-h-1]))
            forecasted_values = mean_values[-h-1:]
            plot_values = mean_values[-h-past_values:]
            plot_index = date_index[-h-past_values:]

            plt.figure(figsize=figsize)
            if intervals == True:
                alpha =[0.15*i/float(100) for i in range(50,12,-2)]
                for count in range(9):
                    plt.fill_between(date_index[-h-1:], error_bars[count], error_bars[-count],
                        alpha=alpha[count])     
            plt.plot(plot_index,plot_values)
            plt.title("Forecast for " + self.data_name)
            plt.xlabel("Time")
            plt.ylabel(self.data_name)
            plt.show()