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 95% prediction intervals for the forecast?

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

def plot_predict(self,h=5,past_values=20,intervals=True,**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 95% prediction intervals for the forecast?

        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:
            # Retrieve data, dates and (transformed) latent variables
            scale, shape, skewness = scale, shape, skewness = self._get_scale_and_shape(self.latent_variables.get_z_values(transformed=True))

            # Get expected values
            forecasted_values = np.zeros(h)

            for value in range(0,h):
                if value == 0:
                    forecasted_values[value] = self.states[0][-1] + self.states[1][-1]
                else:
                    forecasted_values[value] = forecasted_values[value-1] + self.states[1][-1]

            previous_value = self.data[-1]  
            date_index = self.shift_dates(h)
            simulations = 10000
            sim_vector = np.zeros([simulations,h])

            for n in range(0,simulations):  
                rnd_q = np.random.normal(0,np.sqrt(self.latent_variables.get_z_values(transformed=True)[0]),h)
                rnd_q2 = np.random.normal(0,np.sqrt(self.latent_variables.get_z_values(transformed=True)[1]),h)
                exp_0 = np.zeros(h)
                exp_1 = np.zeros(h)

                for value in range(0,h):
                    if value == 0:
                        exp_0[value] = self.states[1][-1] + self.states[0][-1] + rnd_q[value]
                        exp_1[value] = self.states[1][-1] + rnd_q2[value]
                    else:
                        exp_0[value] = exp_0[value-1] + exp_1[value-1] + rnd_q[value]
                        exp_1[value] = exp_1[value-1] + rnd_q2[value]

                sim_vector[n] = self.family.draw_variable(loc=self.link(exp_0),shape=shape,scale=scale,skewness=skewness,nsims=exp_0.shape[0])

            sim_vector = np.transpose(sim_vector)
            forecasted_values = self.link(forecasted_values)

            plt.figure(figsize=figsize) 

            if intervals == True:
                plt.fill_between(date_index[-h-1:], np.insert([np.percentile(i,5) for i in sim_vector],0,previous_value), 
                    np.insert([np.percentile(i,95) for i in sim_vector],0,previous_value), alpha=0.2,label="95 C.I.")   

            plot_values = np.append(self.data[-past_values:],forecasted_values)
            plot_index = date_index[-h-past_values:]

            plt.plot(plot_index,plot_values,label=self.data_name)
            plt.title("Forecast for " + self.data_name)
            plt.xlabel("Time")
            plt.ylabel(self.data_name)
            plt.show()