STUDYING THE EFFECTS OF THE SKEWNESS OF INTER-ARRIVAL AND SERVICE TIMES ON THE PROBABILITY DISTRIBUTION OF WAITING TIMES

Previous studies have shown that the mean queue length of a GI/G/1 system is significantly influenced by the skewness of inter-arrival times, but not by the skewness of service times. These results are limited because all the distributions considered in previous studies were positively skewed. To address this limitation, this paper investigates the effects of the skewness of inter-arrival and service times on the probability distribution of waiting times, when a negatively skewed distribution is used to model inter-arrival and service times. Subsequent to a series of experiments on a GI/G/1 queue using discrete-event simulation, results have shown that the lowest mean waiting time and the lowest variance of waiting times can be attained with a combination of positive inter-arrival skewness and negative service skewness. Results also show an interesting effect of the skewness of service times in the probability of no-delay in environments with a higher utilization factor.


INTRODUCTION
The successful management of waiting time has been important to maintain the competitive advantage of manufacturing companies (de Treville et al., 2004;Tersine & Hummingbird, 1995) as it is a principal component of the manufacturing flow time. Today's business environment has *Corresponding author evolved into a customer-oriented market, where a significant amount of production is made-toorder (Romero-Silva et al., 2016). This change has created a greater need to accurately estimate an order's flow time to carry out various levels of planning throughout the company and reach target service levels.
Queueing theory has been widely used to study waiting times in stochastic environments by modelling manufacturing and service processes as queueing systems. One of the most relevant conclusions of queueing theory literature (Hopp & Spearman, 2000) is that as a system's utilization increases, mean waiting time increases. Furthermore, as either arrival or service time variability increases, mean waiting time also increases (Buzacott & Shanthikumar, 1993).
Other, less frequently mentioned conclusions are that higher-order parameters of inter-arrival and service time distributions, e.g. skewness and kurtosis, influence the mean waiting time of a GI/G/1 queue (Gross & Juttijudata, 1997;Lemoine, 1976), and that the shape of input distributions also influences the resulting mean waiting time (Gross, 1999). Skewness measures how symmetrical the shape of a distribution is. It is very relevant for probability distributions other than the normal distribution, as most production times distributions are skewed (Dudley, 1963;Inman, 1999). A positive skewness indicates that the mass of the distribution is concentrated towards the left of the distribution and that the right tail is longer than the left one; whereas a negative skewness suggests that the mass is concentrated towards the right and that the longer tail is left-located. On the other hand, kurtosis measures the level of mass concentration towards the tails of the distribution.
Regarding skewness, Whitt (1984a) demonstrated that the mean waiting time of a H 2 /H 2 /1 queue (inter-arrival and service times modelled with a second-order hyperexponential distribution (Tarasov, 2016) with single server) is significantly influenced by the skewness of inter-arrival times, whereas the influence of the skewness of the service times was minimal.
Results from Johnson (1993) suggest that having a good approximation for the distribution of inter-arrival times (up to the third moment) is more critical for reducing the error of waiting time estimation than having a good approximation for the distribution of service times, suggesting a higher influence of inter-arrival times on the mean waiting time. However, she suggested that more research is needed regarding the influence of service time skewness on mean queue length approximations.
Despite their value, the results from the studies of Whitt and Johnson are limited since the probability distributions used were all positively skewed and had high coefficient of variation values, which is a characteristic not present in many real production systems (Doerr et al., 2004). Therefore, this paper is concerned with investigating whether the influence of the skewness of service times on the mean waiting time is more significant than previously suggested. In addition, this paper will investigate how the skewness of inter-arrival and service times can further influence the distribution of waiting times, beyond the mean, to significantly impact performance. This study will address the following questions: 1. What is the influence of inter-arrival and service time skewness on the mean waiting time for queues with positively and negatively skewed distributions?
2. How does the skewness of inter-arrival and service times influence the distribution of waiting times?
This topic is investigated through the analysis of the waiting times of a single-server queue, classified as GI/G/1. A GI/G/1 queue entails a single server with general service times and independent, general inter-arrival times distributions. This paper contributes to the body of knowledge by assessing the influence of inter-arrival and service time skewness on the distribution of waiting times. In addition, the managerial implications of this research have the potential to contribute to enhanced production performance by providing additional understanding into how to reduce waiting times.
This paper is organized as follows. Section 2 presents a relevant literature review on the estimation of waiting time parameters and its probability distribution. Section 3 describes the methodology and experimental design used for this research. Section 4 presents the results of the study, while the results discussion is shown in Section 5. Finally, Section 6 presents the conclusions of the study.

LITERATURE REVIEW
Studies investigating waiting time performance have used two main approaches-analytical methods and discrete-event simulation. Most previous literature has applied queueing theory methods to obtain exact solutions to estimate the waiting time probability distributions for PH/PH/1 (inter-arrival and service times modelled as phase-type distributions with single server) queues (Latouche & Ramaswami, 1999a), using matrix analytic methods (Latouche & Ramaswami, 1999b).
In addition, a number of approximations have been obtained for a variety of queueing systems. Buzacott and Shanthikumar (1993) and Bolch et al. (1998) concisely presented a summary of the most relevant GI/G/1 queue approximations. Myskja (1990) presented a list of approximations for the GI/G/1 queue where it can be seen that most approximations consider the coefficient of variation of inter-arrival times (CV A ), the coefficient of variation of service times (CV S ) and the utilization of the server (ρ) as the three most critical factors affecting the mean waiting time (W). Furthermore, Bhat (1993) suggested an approximation for the variance of the waiting time of a GI/G/1 queue where the skewness of the distribution of inter-arrival times (γ A ) and the skewness of the distribution of service times (γ S ) are critical factors.
Based on this previous work, we can see that the k and k+1 moments of the distributions of interarrival and service times are the factors that most influence the kth moment of the mean waiting time, apart from the utilization factor.
Nevertheless, as Lemoine (1976) proved that higher-order moments also have an influence on mean waiting times, a number of papers have since studied the influence of skewness on the mean waiting time. For example, Klincewicz and Whitt (1984) and Johnson and Taafe (1991) showed that the effect of γ A on W is higher as CV A increases, and lower when ρ increases in a GI/M/1 queue. In addition, Sahin and Perrakis (1976) and Whitt (1984a) suggested that the impact of γ A on W of a GI/G/1 queue is significantly higher than the impact of γ S on W, and that the variance of the waiting times is significantly affected by both γ A and γ S .
Lau and Martin (1987) suggested that lower values of processing times' skewness result in higher throughput. They also suggested that the effect of kurtosis on throughput is a complex behaviour to model because the impact of kurtosis depends on its interaction with other factors, such as CV S , γ A and γ S .
Johnson (1993) also investigated the effect that two and three moment approximations have on the estimation of W. She found that, although both γ A and γ S have an effect on W, the highest effect is caused by γ A , and that the effect of γ S is reduced when either CV A or ρ increase. She concluded that more analytical and empirical investigations are needed to extend the conclusions of her study regarding the effect of good service-time approximations.
Powell and Pyke (1994) suggested that negative skewness of processing times results in higher throughput. Furthermore, they show that the effect of kurtosis is difficult to assess. However, they suggest that kurtosis could have an effect on mean throughput up to 5% and that the error of throughput estimation could be as high as 28% when only the first two moments of the distribution of processing times are considered.
More recently, Wu et al. (2018) developed a three moment approximation for the GI/G/1 queue. They stated that a two-moment approximation should be used when CV A is less than one. If CV A is equal or higher than one, they suggested using a three-moment approximation.
In addition, Romero-Silva et al. (2019) studied the effect of skewness on the inter-departure times distribution. They found that combination of negative inter-arrival skewness and positive service skewness reduced the coefficient of variation of inter-departure times, as well as the absolute values for the Lag-1 inter-departures autocorrelation, whereas the opposite was true for a combination of positive inter-arrival skewness and negative service skewness.
In an alternate perspective to the above methodologies, Brandwajn and Begin (2009) studied the effects of the higher order properties of the service times' distribution on the distribution of the number of customers waiting in an M/G/1 queue. They found that the probability distribution tail of the number of customers waiting in queue is influenced by γ S since a higher value of γ S resulted in a shorter tail, reducing the probability of finding a higher number of customers in the system at any moment.
Therefore, the purpose of this paper is to extend the conclusions suggested by previous studies regarding the effects of the skewness of the distributions of inter-arrival and service times on the distribution of waiting times. Specifically, this paper investigates whether the impact of the skewness of service times remains minimal when negatively skewed data with low values of CV are considered.

3 METHODOLOGY
To study whether the skewness of service times has a significant impact on the distribution of waiting times, negatively skewed distributions need to be considered, since previous studies only considered positively skewed distributions. Additionally, as it was important to isolate the effect of skewness on the waiting times, we alternate between positively and negatively skewed random variates, without varying the values of the mean, variance and kurtosis.
Therefore, the recommendation of Khalil et al. (2008) to model inter-arrival and service times with the triangular distribution was followed, as its parameters allow having positive and negative skewness in different experiments, while keeping the values of the mean, variance and kurtosis constant (kurtosis for the triangular distribution equals - 3 5 , irrespective of the distribution parameters). This particularity allows the experimental design to exclusively assess the impact of the skewness on performance, as previous studies have done (see, e.g. Romero-Silva et al., 2019). In addition, the triangular distribution can be used to model simple estimates (Law, 2014) of real random variates as the distributions of some real production processes have low squared coefficient of variation values (Doerr et al., 2004).
Describing the variability of waiting times with exact formulas, whenever non-phase-type distributions are involved, such as the triangular distribution, is challenging (Wolff, 1970). Discrete-Event Simulation (DES) has been used to address this limitation because it is a methodology that can reproduce the behavior of a queueing network better than other techniques (Gue and Kim, 2012; Lagershausen & Tan, 2015). Hence, DES was used as the modelling paradigm for this experimental study.
Four waiting time responses were selected in this study to capture the effect of γ A and γ S on the distribution of waiting times. Typical moment-related measures were chosen as responses, namely, the mean (W), variance (σ 2 W ) and skewness (γ W ) of the distribution of waiting times, as well as the probability that the customer does not wait in queue (PW0). However, the actual responses W, σ 2 W , γ W and PW0 reported in this study are the averages of the values of each moment-related measure throughout the 200 replications. In Table 1, we can find a reference  table of all the variables  To determine the length of the transient, Welch's (1983) method was used. Five initial replications were carried out to measure the number of system customers, leading to the conclusion that around 1,000 initial customers were needed to reach the steady-state. This conclusion was reached by comparing the expected number of customers for the given experimental factors, which were estimated by using Hopp and Spearman's (2000) approximation formula for mean number of customers in a GI/G/1 steady-state queue, against the moving averages obtained in Thus, if w kij is the waiting time of customer j in replication i, and experiment k; N is the number of customers per replication (10,000), and R is the number of replications (200); then the mean waiting time of all customers in replication i for experiment k (ẃ ki ) was calculated as: w kij for all i = 1, 2, . . . , R

7
In addition, the mean waiting time of all replications in experiment k (W k ), the waiting time variance of all replications in experiment k (σ 2 wk ), and the waiting time skewness of all replications in experiment k (γ Wk ) were calculated in the following manner: ki , for all k = 1, 2, . . . , 252 Furthermore, the probability that a customer does not wait for service (no-delay) in replication i, experiment k, (p0 ki ) was calculated as follows: p0 ki = number of customers with zero waiting time ∈ the ith replication N Thus, the probability of no-delay for all customers in experiment k was estimated as: Four factors were considered as independent variables: γ A , γ S , CV of the system and ρ , as CV A , CV S and ρ have been shown to have an interaction with the effects of γ A and γ S . Note that since the objective is to identify the interactions of the probability distributions and not the interactions among different levels of general variability, CV A was set equal to CV S , so that the results are not influenced by variability interactions.
Three factor levels were defined for γ A and γ S , representing a negative, zero and positive skewness, while two different values of the system's CV were defined to further study the interaction between skewness, CV and ρ. Moreover, several values of ρ were considered to represent low, medium and very high traffic intensities. Thus, a total of 252 different experiments were run.   distribution were: mean = 3.67, standard deviation = 0.85, and skewness = 0.42. In order to consider an equivalent triangular distribution with equal mean and standard deviation but negative skewness, we calculated the minimum value of parameter a that will result in a mean of 3.67 and a standard deviation of 0.85, by using a second order equation that equalled the coefficient of variation of the positively and negatively skewed triangular distributions. With the resulting value (a = 1.26), the maximum value of parameter m that will result in a mean of 3.67 and a standard deviation of 0.85 was calculated in order to attain a negative skewness. Finally, as the mean and standard deviation needed to be equal, there was only one possible value of b, considering the previous results obtained from a and m. These calculations resulted in a triangular distribution with a skewness of -0.57.
To model a non-skewed triangular distribution for inter-arrival times, we considered that the mode was equal to the mean of the triangular distribution (m = 3.67) and, based on that initial value, we calculated the value of b that would result in a mean of 3.67 and a standard deviation of 0.85. Then we calculated the only possible value of a, given the values of m and b. A graphical representation of the resulting triangular density probability functions is presented in Figure 1 (a).  In addition, based on the original parameters of Khalil et al. (2008) of a negatively skewed triangular distribution (a = 0, m = 3, b = 4), the resulting standardized moments of the triangular distribution were: mean = 2.33, standard deviation = 0.85, and skewness = -0.42. Taking into consideration that the mean and standard deviation should remain constant, the parameters for a positively and non-skewed distribution were calculated in a similar manner as before. A graphical representation of the resulting triangular density functions resulting from this exercise can be found in Figure 1 (b).
On the other hand, the parameters for the triangular distributions modelling the service times were calculated by taking into account the fact that for a single server queue, ρ = mean service time/mean inter-arrival time. Therefore, the mean service time was calculated based on the experimental values of ρ and the corresponding experimental mean inter-arrival time. With the mean service time value and, given that CV S should be equal to CV A (by experimental design), we calculated the standard deviation for the service times. As the mean and the standard deviation of service times were now known, the parameters of the triangular distribution (a, m and b) were then calculated with the same procedure used to calculate the parameters of the distributions of inter-arrival times for positively, negatively and non-skewed distributions. The full experimental values of inter-arrival and service times stemming from this procedure are shown in the Appendix, Table A1.
Groups of nine experiments with the same values of CV and ρ, but different values of γ A and γ S , were constituted to directly assess the impact of γ A and γ S on the distribution of waiting times, resulting in a total of 28 groups. For instance, Group A was constituted of experiments where CV = 0.2318 and ρ = 0.1 (see Table A1 in the Appendix).
ANOVA tests for each group were run to gauge whether a difference in the values of γ A and γ S created statistically significant differences between the responses. In addition, Duncan's multiple range tests (Duncan, 1955) were carried out for each group to see whether statistical differences between experiments within the group existed.
Furthermore, for each experimental response, an intra-group variation (IGV) calculation was carried out by dividing the actual value of the response in each experiment by the average of the results of the nine experiments included in the corresponding group, resulting in five additional responses (IGV W , IGV σ 2 , IGV γ and IGV PW0 ), whose values oscillate around 1. As IGV values represent a deviation from the average of the group, they can be used to compare inter-group experiments.
For example, to calculate the intra group variation for the mean waiting time of experiment k (IGV Wk ), included in Group x, the following formula was used: Therefore, IGV calculations were used to estimate the statistical impact of each independent variable, as well as the interactions between the variables on each of the five responses, by carrying out ANOVA tests of the impact of single factors and their interactions. All statistical tests were run using R (The R Foundation, 2019).
Finally, to simplify the analysis of the results, correlation tests between pairs of responses were completed, resulting in high and statistically significant correlations between W and σ 2 W ; therefore, only the results for W, γ W and PW0 were further analyzed.

RESULTS
After conducting the ANOVA and Duncan's multiple range tests, some groups of experiments with ρ > 0.90 yielded non-statistically significant differences among experiments while some experiments with ρ < 0.50 resulted in a W = 0. Therefore, the results presented in this section do not include some of these groups of experiments. Experimental results for these experiments can be found in Table A3 in the Appendix, whereas Table A4 in the Appendix shows the Duncan's test results.

Mean waiting time
The impact of γ A and γ S on W at different levels of ρ is shown in Figure 2 and Figure 3. From those figures it can be seen that the effect of the skewness of inter-arrival and service times on the mean waiting time diminishes as the server's utilization increases because the relative differences among the experiments are reduced with increasing values of server's utilization. The results shown in Figure 2 and Figure 3 support previous conclusions from other studies (Johnson, 1993; Whitt, 1984a) as they demonstrate the impact of γ A on W, i.e. positive γ A reduces W. However, they also show that γ S has a significant influence on W, as the differences between the values of W in experiments with different γ S but equal γ A are considerable. For example, for experiments with CV = 0.2318, ρ = 0.8 and positive γ A (see Figure 2), the resulting W values are 0.2959, 0.3115 and 0.3291 for negative, zero and positive γ S values (see Table A3 in the Appendix), respectively. This entails more than 10% deviation between the maximum and minimum values, when compared with the average of these three results. Moreover, the 5% difference between the experiment with positive γ S and that with zero γ S is also sizeable, showing that the impact of γ S on W is significant even for moderate changes. The results also show that for every experimental value of γ A, a negative γ S produces the best W values.
The effect of both γ A and γ S on W seems to decrease with a higher CV value, as the values of IGV W for experiments with a CV = 0.3642 suggest (see Figure 3). Figure 3 shows how lower γ S values decrease W. Furthermore, the results presented in Figure 2 and Figure 3 indicate that a combination of positive γ A and negative γ S produces the best results in terms of W.
The ANOVA test for the effects of the experimental factors on IGV W confirms the aforementioned impact of γ A on W, as this factor, together with its interactions with ρ or CV, has the highest percentage contribution to the variation of IGV W , as shown in Table 3. It can also be observed that the effect of γ S on IGV W is significant, albeit not as important as that of γ A .    Figure 4 exhibit the results of IGV γ (intra-group variation of waiting time skewness). The top-left plot indicates that a positive γ A causes an increase in γ W , which is directly related to W results, as a positive γ A decreases the mean waiting times, and a higher CV value diminishes the impact of γ A on both γ W and W. This is an opposite result from the effect of the interaction between CV and γ S on γ W (top, middle-left plot), where a higher CV value creates a higher effect of γ S on γ W . Interestingly, the effect of γ S on γ W is different for experiments with different CV values. Positive γ S decreased γ W in experiments with CV = 0.2318, whereas a positive γ S increased γ W in experiments with CV = 0.3642.

Waiting time skewness and probability of no-delay
The plot (top, middle-right of Figure 4) of the interaction between ρ and γ A suggests that with higher values of ρ, the effect of γ A on γ W drastically diminishes. More interestingly, the inter-  action between ρ and γ S (top-right plot in Figure 4) shows that the effect of γ S on γ W stabilizes and reverses as ρ increases, at least for ρ ≤ 0.85.
In addition, the bottom row of Figure 4 depicts the effect of various factor interactions on IGV PW0 (intra-group variation of no-delay probability) results. These results should be the opposite as the results from W, as higher W values will imply lower PW0. Nevertheless, various results regarding the effect of the factors on PW0 suggest a different behavior. For instance, a positive γ A , which reduces W, also diminishes PW0 with a low CV value; a negative γ S , which reduces W, also diminishes PW0 for the two CV values considered; a positive γ A reduces PW0 for 0.70 < ρ < 0.95 and negative γ A increases PW0 for the same range of ρ values; and negative γ S reduces PW0 for ρ > 0.70 while positive γ S increases PW0 for the same range.
Furthermore, the interactions of γ A and γ S with ρ on IGV PW0 show very intriguing behavior as the curves of the interaction between γ A and ρ suggest a non-linear behavior: the effect of zero γ A shows a decreasing pattern with increasing ρ values, and the effect of γ S on PW0 seems to increase as ρ increases, which contrasts with the results concerned with W, where ρ created a limit on the effect of the skewness on W.
In fact, the effect of γ S on PW0 appears to be very significant, as the percentage contribution of γ S (including its interaction with ρ) to the variation of IGV PW0 is very high (more than 64%, as is shown in Table 4). Interestingly, the ANOVA test for the effects on IGV PW0 also suggests a relevant interaction between γ A and γ S on the values of PW0, an interaction that shows that γ S behaves as expected with positive γ A , i.e. negative γ S increases P W0 as it reduces W, but for negative or zero γ A , a negative γ S reduces PW0 (see Figure A1 in the Appendix). To get a comprehensive understanding of the effect of skewness on waiting times, we decided to further investigate the effect of γ A and γ S on the curve of the distribution of non-zero waiting times. Figure 5 gives an example of these effects by showing three plots comparing combinations of contrasting values of γ A and γ S for different levels of ρ. These plots confirm the main findings of this study, namely, a combination of positive γ A and negative γ S (dotted line) reduces W by shortening the tail of the distribution of waiting times and increasing γ W and PW0. Contrarily, a combination of negative γ A and positive γ S (solid line) increases W by extending the tail of the distributions of waiting times and reducing γ W and PW0. However, this difference is significantly reduced with increasing values of ρ (right plot).

DISCUSSION
Previous studies have suggested that the mean waiting time is minimally affected by the skewness of the distribution of service times (Johnson, 1993;Whitt, 1984a). Adding nuance to this understanding, results obtained from this study show that, although the impact of the skewness of inter-arrival times is higher than that of the skewness of service times, the skewness of service times should also be considered when estimating the mean waiting time, particularly when opposing values of inter-arrival and service skewness are present, as the lowest values for the mean waiting time, for all groups of experiments, were found in experiments with positive values of Thus, to answer the first research question of this study (What is the influence of the skewness of inter-arrival and service times on the mean waiting time for queues with positively and negatively skewed distributions?), it was found that the effect of inter-arrival skewness on the mean waiting time is negative, i.e. the higher the value of γ A , the lower the value of the mean waiting time.
Contrarily, the effect of service time skewness is positive: the higher the value of γ S , the higher the value of the mean waiting time. However, both of these effects are limited when the server's utilization increases as statistically significant differences between experiments with varying values of γ A and γ S were not found for scenarios with a higher system load (ρ > 0.90). This a result that is consistent with heavy-traffic behavior since higher moments than the first two have little effect on the performance of queueing systems when the utilization factor is high ( To answer the second research question of this study (How does the skewness of inter-arrival and service times influence the characteristics of the distribution of waiting times?), it was found that a combination of positive γ A and negative γ S also reduces the variance of waiting times, while increasing the skewness of the distribution of waiting times, resulting in a "taller" curve with a shorter tail. The combination of negative γ A and positive γ S resulted in the opposite results. Again, these results are highly dependent on the utilization factor, as higher values of ρ resulted in opposite effects of γ A and γ S regarding γ W and PW0.
It is interesting to note that the experimental results showed that even though γ A is more significant in reducing the mean waiting times of a GI/G/1 queue than γ S , the effect of γ S on the resulting PW0 is significantly bigger than that of γ A , indicating that practitioners should take into account both γ A and γ S to increase different performance measurements related to waiting times. Moreover, we found a counterintuitive effect of γ S on PW0 since PW0 was reduced with a negative γ S , contradicting expected results that PW0 should increase with negative γ S as a negative γ S reduces W for the same CV value. This result can be compared with another counterintuitive result presented by Whitt (1984b) where it was shown that PW0 was reduced when CV S decreased and CV A > 1 in GI/G/1 queues. Thus, as Whitt explains, more variability caused by increasing values of γ S , instead of increasing values of CV S , produces higher variance in the resulting waiting times with the associated higher mean waiting times, which in turn results in longer tails of the distribution of waiting times and a higher probability of no-delay (the left tail of the distribution) when ρ > 0.70. Moreover, the effect of γ S on IGV PW0 was found to be more important with increasing ρ values, a result that contrasts with the limit that ρ imposes on the effect that γ A and γ S have on W. Whitt (1984a) reported some differing conclusions on the effect of γ S on mean queue length, as he showed that for a CV 2 A value equal to 2, a higher value of γ S created a lower mean queue length, although the variation was very small. However, when carrying out a series of calculations for a system with Erlang inter-arrivals and a CV 2 A equal to 0.5, he found that higher values of γ S created higher values of mean queue length, a result that corresponds to the results of the current study.
These divergent research outcomes suggest that more studies are needed to discover whether contrasting values of CV A have an influence on the effect of γ S on the mean queue length and the mean waiting time, as it would seem that when CV A is lower than 1, higher γ S results in higher mean waiting times; whereas when CV A is higher than 1, higher γ S results in lower mean waiting times.
On the other hand, other study results here differ from those provided by Brandwajn and Begin (2009) regarding the tail of the distribution of the M/G/1 queue length, as they found that a higher value of γ S caused a shorter tail when a CV higher than 1 was considered. Experimental results from the current study show the opposite effect, as higher values of γ S cause a larger tail for the waiting times distribution. Thus, results from the current study suggest that varying the skewness of input distributions can be a good management tool to decrease the risk of having a customer wait a very long time or not at all.
Furthermore, Whitt's results (1984a) showed that the maximum percentage difference in mean waiting times that can be attained by varying positive values of γ A in a H 2 /H 2 /1 queue with ρ = 0.70 reached 30% improvement. The results of the current study suggest that by varying both γ A and γ S the biggest difference in mean waiting times between experiments reached 57% when CV = 0.2318 and ρ = 0.70, whereas when only γ S was modified, the improvement approached 19% in experiments with positive γ A . This shows the importance of considering both γ A and γ S and the impact of γ S on the performance of single server queues with triangular distributions.
Moving forward, results from this paper suggest that managers and researchers of GI/G/1 queues with low CV A and CV S values should consider the effect of service time skewness on the estimation of the mean waiting time, the probability of no-delay and waiting times distribution tail in order to get accurate estimations. Therefore, approximation methods estimating mean waiting times for GI/G/1 queues with CV < 1 should consider both inter-arrival and service skewness in their approximation procedures (see, e.g. Whitt, 1989; Wu et al., 2018).

Limitations of the study and future research
Although the selection of DES for this study was well-suited for this application, as with all modelling methodologies, it contains assumptions and limitations which detract from its true replication of reality. Moreover, experimental analysis, although enhancing control and providing structure, limits the generalizability of the conclusions provided by the experiments, in comparison with analytical approaches.
Given these parameters, the results of this study are limited to only the particular model studied here, although the results do provide valuable practical and theoretical insights on the effects of the skewness of inter-arrival and service times on the probability distribution of waiting times. Future research streams are needed to determine if the insights found here are complementary and applicable for multiple server (GI/G/m) systems, serial lines and single-server queues with high values of CV A and CV S .
In addition, this study is based on a limited sample of experiments, with few variations among the values of the factors, and only investigated low-variability systems with triangular distributions, by incorporating only two CV values. Thus, more research that considers probability distributions with contrasting values of CV is needed, as well as additional studies that consider much higher and lower γ A and γ S values in order to study the impact that different magnitudes of skewness have on the performance of single-server queues. However, since this is the first study that explored the effect of negatively skewed distributions on mean waiting times, it provides us with further motivation to investigate a wider range of skewness values, as it was found that varying the service time skewness could lead to an improvement of up to 19% in the performance of the queue. Furthermore, these results remain applicable for some realistic contexts in terms of CV values, as it has been shown that some production processes have low values of CV (Doerr et al., 2004;Inman, 1999) and very similar skewness magnitudes (Slack, 1982) than the ones considered in this paper.
More studies are needed to investigate the actual impact of the kurtosis on the performance of single queues since previous studies (Powell & Pyke, 1994) have shown that it could impact mean throughput up to 5% in saturated production lines and since the experimental design of this current study specifically omitted this impact by using the triangular distribution, which has constant kurtosis, with the intention of exclusively assessing the impact of skewness.
Lastly, it might be argued that the applicability of the results is limited because the combination of inter-arrival and service skewness that reduced the mean waiting times in this study can be difficult to find in real contexts, as negatively skewed service times distributions are rare. However, this study shows that the difference in mean waiting times between an experiment with positive service skewness and zero service skewness, both commonly present in real scenarios (Inman, 1999), is significant and it also reduces mean waiting time. Furthermore, it has been shown (Dudley, 1963) that changing from positively skewed service processes to a process with no skewness can be attained in practice.

CONCLUSIONS
This paper investigated the interactions among the probability distributions of inter-arrival and service times and their effects on the distribution of waiting times, and contributes to the current literature by filling a long-standing gap in previous research. Experimental analysis of a singleserver queue with triangular distributions of inter-arrival and service times showed that the mean waiting time was influenced by both inter-arrival and service skewness, and that such influence is limited by increasing CV and server's utilization values.
Results suggest that a combination of positive inter-arrival skewness and negative service skewness result in smaller values of mean waiting times, a reduced variance of waiting times, a higher probability of no-delay, and an incremental skewness and shorter tail for the distribution of waiting times, increasing the performance of the GI/G/1 queue when the CV of inter-arrival times is lower than 1 and the server's utilization is low or moderate. In contrast, negative service skewness in higher traffic environments will result in a lower probability of no-delay while at the same time producing lower mean waiting times and variance of waiting times.
Future research is needed to more comprehensively study how the coefficient of variation of inter-arrival times influences the effect service skewness on the mean waiting time, as previous results considering a CV of inter-arrival times higher than 1 have found contrary results to the ones reported in the current study.  Table A1 -Experimental values for the study's experiments (see Table A2 for the variables' abbreviations).      Figure A1 -Plot of the effect of the interaction between γ A and γ S on PW0.