code: def first_order_score(y, mean, scale, shape, skewness):
          
          m1 = (np.sqrt(shape)*sp.gamma((shape-1.0)/2.0))/(np.sqrt(np.pi)*sp.gamma(shape/2.0))
          mean = mean + (skewness - (1.0/skewness))*scale*m1
          if (y-mean)>=0:
              return ((shape+1)/shape)*(y-mean)/(np.power(skewness*scale,2) + (np.power(y-mean,2)/shape))
          else:
              return ((shape+1)/shape)*(y-mean)/(np.power(scale,2) + (np.power(skewness*(y-mean),2)/shape))
docs: GAS Skew t Update term using gradient only - native Python function
  
          Parameters
          ----------
          y : float
              datapoint for the time series
  
          mean : float
              location parameter for the Skew t distribution
  
          scale : float
              scale parameter for the Skew t distribution
  
          shape : float
              tail thickness parameter for the Skew t distribution
  
          skewness : float
              skewness parameter for the Skew t distribution
  
          Returns
          ----------
          - Score of the Skew t family