Again

BookTwo <- read.csv("FiveYears.csv", header = TRUE, stringsAsFactors = FALSE)
kable(head(BookTwo))

Item_Code	Description	X1.1.2012	X2.1.2012	X3.1.2012	X4.1.2012	X5.1.2012	X6.1.2012	X7.1.2012	X8.1.2012	X9.1.2012	X10.1.2012	X11.1.2012	X12.1.2012	X1.1.2013	X2.1.2013	X3.1.2013	X4.1.2013	X5.1.2013	X6.1.2013	X7.1.2013	X8.1.2013	X9.1.2013	X10.1.2013	X11.1.2013	X12.1.2013	X1.1.2014	X2.1.2014	X3.1.2014	X4.1.2014	X5.1.2014	X6.1.2014	X7.1.2014	X8.1.2014	X9.1.2014	X10.1.2014	X11.1.2014	X12.1.2014	X1.1.2015	X2.1.2015	X3.1.2015	X4.1.2015	X5.1.2015	X6.1.2015	X7.1.2015	X8.1.2015	X9.1.2015	X10.1.2015	X11.1.2015	X12.1.2015	X1.12.2016	X2.12.2016	X3.12.2016	X4.12.2016	X5.12.2016	X6.12.2016	X7.12.2016	X8.12.2016	X9.12.2016	X10.12.2016
SS02CLR3660	3/32 CLEAR 36X60	18	3	5	1	1	2	1	6	0	2	5	6	12	8	3	7	3	2	0	108	1	3	0	3	2	3	1	0	1	0	0	0	0	0	17	0	0	0	0	0	0	0	0	0	0	0	1	1	0	0	0	1	0	0	0	0	1	0
SS02CLR4884	3/32 CLEAR 48X84	2	0	4	6	2	7	0	1	0	92	0	1	3	5	1	0	0	1	1	35	157	1	0	1	1	0	0	1	0	1	0	2	7	0	0	2	0	0	0	0	213	0	2	1	0	0	1	0	1	4	7	4	152	2	0	0	1	1
SS02CLR7284	3/32 CLEAR 72 X 84	468	538	687	932	1,162	1,067	1,131	1,201	1,183	1,082	793	505	498	382	401	752	955	878	999	1,053	1,007	1,283	753	396	397	381	422	742	1,028	1,168	1,216	1,179	1,329	1,161	575	622	347	296	566	820	919	1,379	1,159	1,132	1,119	1,238	904	699	512	471	744	789	943	1,341	973	1,255	1,146	1,045
SS02LOE7284	3/32 ENEG ADV LOW-E 72X84	0	0	0	0	64	186	176	185	186	191	141	72	92	62	64	139	208	162	188	198	194	217	124	89	71	84	53	140	203	214	221	226	262	239	115	100	64	45	97	195	192	257	235	250	275	256	208	134	102	82	173	146	232	263	215	262	261	248
SS03ACC7284	1/8 ACCLIMATE 72X84	30	33	44	48	117	93	133	126	140	101	48	37	29	24	23	40	79	39	57	28	18	15	9	6	17	5	6	8	15	20	21	9	8	19	14	20	19	8	27	15	44	26	20	18	20	15	14	10	4	6	8	3	11	15	8	24	9	7
SS03BRN7284	1/8 BRONZE 72X84	0	0	26	40	0	0	0	0	0	0	0	0	0	0	0	0	0	0	23	42	28	24	6	0	0	0	0	0	0	0	0	0	30	23	7	0	0	70	226	23	20	37	32	65	42	36	46	62	59	13	35	59	32	0	0	0	0	0

for(i in c(3:20:ncol(BookTwo))){
  BookTwo[,i] <- as.numeric(BookTwo[,i])}
BookTwo[is.na(BookTwo)] <- 0 
    BookTwo2 <- gather(BookTwo,Date, Inventory_Used, -c(Item_Code, Description))

   
    BookTwo2$Date <- c("X1.1.2012", "X2.1.2012", "X3.1.2012", "X4.1.2012", "X5.1.2012", "X6.1.2012", "X7.1.2012", "X8.1.2012", "X9.1.2012", "X10.1.2012", "X11.1.2012", "X12.1.2012","X1.1.2013", "X2.1.2013", "X3.1.2013", "X4.1.2013", "X5.1.2013", "X6.1.2013", "X7.1.2013", "X8.1.2013", "X9.1.2013", "X10.1.2013", "X11.1.2013", "X12.1.2013","X1.1.2014", "X2.1.2014", "X3.1.2014", "X4.1.2014", "X5.1.2014", "X6.1.2014", "X7.1.2014", "X8.1.2014", "X9.1.2014", "X10.1.2014", "X11.1.2014", "X12.1.2014","X1.1.2015", "X2.1.2015", "X3.1.2015", "X4.1.2015", "X5.1.2015", "X6.1.2015", "X7.1.2015", "X8.1.2015", "X9.1.2015", "X10.1.2015", "X11.1.2015", "X12.1.2015","X1.1.2016", "X2.1.2016", "X3.1.2016", "X4.1.2016", "X5.1.2016", "X6.1.2016", "X7.1.2016", "X8.1.2016", "X9.1.2016", "X10.1.2016")
    
    BookTwo2$Date <- as.Date(substring(BookTwo2$Date, 2), "%m.%d.%Y")
    
    kable(head(BookTwo2))

Item_Code	Description	Date	Inventory_Used
SS02CLR3660	3/32 CLEAR 36X60	2012-01-01	18
SS02CLR4884	3/32 CLEAR 48X84	2012-02-01	2
SS02CLR7284	3/32 CLEAR 72 X 84	2012-03-01	468
SS02LOE7284	3/32 ENEG ADV LOW-E 72X84	2012-04-01	0
SS03ACC7284	1/8 ACCLIMATE 72X84	2012-05-01	30
SS03BRN7284	1/8 BRONZE 72X84	2012-06-01	0

    str(BookTwo2)

## 'data.frame':    6206 obs. of  4 variables:
##  $ Item_Code     : chr  "SS02CLR3660" "SS02CLR4884" "SS02CLR7284" "SS02LOE7284" ...
##  $ Description   : chr  "3/32 CLEAR 36X60" "3/32 CLEAR 48X84" "3/32 CLEAR 72 X 84" "3/32  ENEG ADV LOW-E 72X84" ...
##  $ Date          : Date, format: "2012-01-01" "2012-02-01" ...
##  $ Inventory_Used: num  18 2 468 0 30 0 3 0 7 0 ...

meanItems <- BookTwo2 %>% select(Item_Code, Description, Inventory_Used) %>%  group_by(Item_Code) %>% summarise(Monthly_Average = mean(Inventory_Used), StDev_Monthly = sd(Inventory_Used), CV = StDev_Monthly/Monthly_Average, Av_Demand_during_L = Monthly_Average*.214, Safety_Stock = StDev_Monthly*2.05*sqrt(.214), Reorder_Level = Av_Demand_during_L+2.05*StDev_Monthly*sqrt(.214), Order_Quantity = sqrt(2*Monthly_Average)) %>% arrange(desc(Monthly_Average)) %>% slice(1:30)

Problem Definition

In production and inventory control, increased accuracy is likely to lead to lower safety stocks. Here the manager and forecaster must weigh the cost of a more sophisticated and more expensive technique against potential savings in inventory costs.Harvard Business Review

An inventory system that prevents stockouts.

*Safety stock

*Forecasting

Decide on a horizon. A month?

*Prediction probability? Looking to prepare for a 20% increase in sales.

*This is time series forecast, we are not looking into factors which may affect the SKU behavior.

Data Gathering

How did we decide what data to pull?

*Join 2015 and 2016 monthly inventory data

*This narrowed our data ser to 175 SKUs

*Calculate the Monthly_Average monthly use of each SKU

*Chose the top 30 SKUs with the highest Monthly_Average

Preliminary Analysis

Graphing the data

Are there consistent patterns?

Is there a significant trend?

Is seasonality important?

Is there evidence of the presence of business cycles? Are there any outliers in the data that need to be explained by those with expert knowledge?

How strong are the relationships among the variables available for analysis?

Find the highest volume SKUs and do inventory management calculations based on historical data.

Equations based on the following assumptions:

1. Continuous Review Policy

2. consistent lead time of 6 days and a service level of 98%.

3. Service level of 98%.

Average Monthly Demand:\[\bar{x}_d\]

Standard Deviation of Monthly Demand:\[\sigma_d\]

Lead Time:\[L = .214\]

Service Level:\[z = 2.05\]

Carrying Cost (negligible):\[k = 1\]

Holding Cost (negligible):\[h = 1\]

Average demand during lead time: \[\bar{x}_L = \bar{x}_d * L\]

Safety Stock:\[z*\sigma_d*\sqrt{L}\]

Reorder Level:\[\bar{x}_L+z*\sigma_d*\sqrt{L}\]

Order Quantity:\[Q = \sqrt\frac{2k*\bar{x}_d}{h}\ \]

Which is the best way to show historical data?

TwoYr222 (660,11) All calculations, all months, top 30

kable(head(meanItems), digits = 2)

Item_Code	Monthly_Average	StDev_Monthly	CV	Av_Demand_during_L	Safety_Stock	Reorder_Level	Order_Quantity
SS05CLR96130	403.69	121.05	0.30	86.39	114.79	201.18	28.41
SS02CLR7284	380.33	361.88	0.95	81.39	343.18	424.57	27.58
SS10CLR96130	327.03	110.00	0.34	69.99	104.31	174.30	25.57
SS06S6T100144	321.69	210.00	0.65	68.84	199.15	267.99	25.36
SS12CLR102130	309.98	138.67	0.45	66.34	131.51	197.84	24.90
SS04CLR96130	225.41	68.36	0.30	48.24	64.83	113.07	21.23

Item_Code	Description	Date	Inventory_Used	Monthly_Average	StDev_Monthly	CV	Av_Demand_during_L	Safety_Stock	Reorder_Level	Order_Quantity	Z_Score
SS02CLR7284	3/32 CLEAR 72 X 84	2012-03-01	468	380.33	361.88	0.95	81.39	343.18	424.57	27.58	0.24
SS02LOE7284	3/32 ENEG ADV LOW-E 72X84	2012-04-01	0	156.17	79.04	0.51	33.42	74.95	108.38	17.67	-1.98
SS03ACC7284	1/8 ACCLIMATE 72X84	2012-05-01	30	31.21	33.66	1.08	6.68	31.92	38.60	7.90	-0.04
SS03CLR7284	1/8 CLEAR 72X84	2012-10-01	0	133.53	312.50	2.34	28.58	296.36	324.93	16.34	-0.43
SS03S6A7284	1/8 SOLARBAN60 72X84	2013-05-01	3	63.71	51.48	0.81	13.63	48.82	62.45	11.29	-1.18
SS03S6T7284	1/8 SOLARBAN60 TC 72X84	2013-06-01	2	61.59	59.73	0.97	13.18	56.64	69.82	11.10	-1.00

I22 <- TwoYr222 %>% filter(Item_Code == "SS06CLR96130") %>% select(Date, Item_Code, Inventory_Used, Z_Score)
XV22<- TwoYr222 %>% filter(Item_Code == "SS03CLR7284") %>% select(Date, Item_Code, Inventory_Used, Z_Score)
XI22<- TwoYr222 %>% filter(Item_Code == "SS02CLR7284") %>% select(Date, Item_Code, Inventory_Used, Z_Score)
IV22<- TwoYr222 %>% filter(Item_Code == "SS06S6T100144") %>% select(Date, Item_Code, Inventory_Used, Z_Score)
V22<- TwoYr222 %>% filter(Item_Code == "SS05CLR96130") %>% select(Date, Item_Code, Inventory_Used, Z_Score)
VI22<- TwoYr222 %>% filter(Item_Code == "SS12CLR102130") %>% select(Date, Item_Code, Inventory_Used, Z_Score)
VII22<- TwoYr222 %>% filter(Item_Code == "SS10CLR96130") %>% select(Date, Item_Code, Inventory_Used, Z_Score)
VIII22<- TwoYr222 %>% filter(Item_Code == "SS03CLR96130") %>% select(Date, Item_Code, Inventory_Used, Z_Score)
IX22<- TwoYr222 %>% filter(Item_Code == "SS06LCL4080") %>% select(Date, Item_Code, Inventory_Used, Z_Score)
X22<- TwoYr222 %>% filter(Item_Code == "SS06STC96130") %>% select(Date, Item_Code, Inventory_Used, Z_Score)
III22<- TwoYr222 %>% filter(Item_Code == "SS06S6A100144") %>% select(Date, Item_Code, Inventory_Used, Z_Score)
XII22<- TwoYr222 %>% filter(Item_Code == "SS03LOE7284") %>% select(Date, Item_Code, Inventory_Used, Z_Score)
XIII22<- TwoYr222 %>% filter(Item_Code == "SS10CLV96130") %>% select(Date, Item_Code, Inventory_Used, Z_Score)
XIV22<- TwoYr222 %>% filter(Item_Code == "SS04CLR96130") %>% select(Date, Item_Code, Inventory_Used, Z_Score)
II22<- TwoYr222 %>% filter(Item_Code == "SS06MCL96130") %>% select(Date, Item_Code, Inventory_Used, Z_Score)
XVI22<- TwoYr222 %>% filter(Item_Code == "SS06LOE96130") %>% select(Date, Item_Code, Inventory_Used, Z_Score)
XVII22<- TwoYr222 %>% filter(Item_Code == "SS02LOE7284") %>% select(Date, Item_Code, Inventory_Used, Z_Score)
XVIII22<- TwoYr222 %>% filter(Item_Code == "SS10SHG96130") %>% select(Date, Item_Code, Inventory_Used, Z_Score)
XIX22<- TwoYr222 %>% filter(Item_Code == "SS03S6T7284") %>% select(Date, Item_Code, Inventory_Used, Z_Score)
XX22<- TwoYr222 %>% filter(Item_Code == "SS03S6A7284") %>% select(Date, Item_Code, Inventory_Used, Z_Score)
XXI22<- TwoYr222 %>% filter(Item_Code == "SS06STP96130") %>% select(Date, Item_Code, Inventory_Used, Z_Score)
XXII22<- TwoYr222 %>% filter(Item_Code == "SS06S7T100144") %>% select(Date, Item_Code, Inventory_Used, Z_Score)
XXIII22<- TwoYr222 %>% filter(Item_Code == "SS05LOE96130") %>% select(Date, Item_Code, Inventory_Used, Z_Score)
XXIV22<- TwoYr222 %>% filter(Item_Code == "SS06GRY96130") %>% select(Date, Item_Code, Inventory_Used, Z_Score)
XXV22<- TwoYr222 %>% filter(Item_Code == "SS03AHT7284") %>% select(Date, Item_Code, Inventory_Used, Z_Score)
XXVI22<- TwoYr222 %>% filter(Item_Code == "SS06BRN96130") %>% select(Date, Item_Code, Inventory_Used, Z_Score)
XXVII22<- TwoYr222 %>% filter(Item_Code == "SS05S6T96130") %>% select(Date, Item_Code, Inventory_Used, Z_Score)
XXVIII22<- TwoYr222 %>% filter(Item_Code == "SS04LOE96130") %>% select(Date, Item_Code, Inventory_Used, Z_Score)
XXIX22<- TwoYr222 %>% filter(Item_Code == "SS10STP96130") %>% select(Date, Item_Code, Inventory_Used, Z_Score)
XXX22<- TwoYr222 %>% filter(Item_Code == "SS12STP96130") %>% select(Date, Item_Code, Inventory_Used, Z_Score)

KeyItems22 <- rbind(I22,II22,III22,IV22,V22,VI22,VII22,VIII22,IX22,X22,XI22,XII22,XIII22,XIV22,XV22,XVI22,XVII22,XVIII22,XIX22,XX22,XXI22,XXII22,XXIII22,XXIV22, XXV22,XXVI22,XXVII22,XXVIII22,XXIX22,XXX22)

ggplot(KeyItems22, aes(x = Date, y = Z_Score)) + geom_line() + facet_wrap(~Item_Code)

ggplot(KeyItems22, aes(x = Date, y = Z_Score)) + geom_smooth() + facet_wrap(~Item_Code)

1/4 CLEAR 96X130

A <- ts(c(1074,1381,1181,1143,1220,1176,1379,1667,1413,1830,1400,1114,1217,1203,1324,1805,1474,1370,2239,1817,1884,2188,2118,2026,1626,1526,1099,1422,1979,1752,1909,2196,2275,2613,2327,2649,2008,1520,1695,2076,1725,2184,2591,2627,2590,2975,2349,2800,1926,1870,2327,2273,2418,2596,2578,3146,2656,2386
), 
        start = c(2012, 1), frequency = 12)

plot(A)

seasonplot(A,ylab="Demand", xlab="Month",
           main="Seasonal plot: 1/4 CLEAR 96X130 Demand",
           year.labels=TRUE, year.labels.left=TRUE, col=1:20, pch=19)

monthplot(A,ylab="Demand",xlab="Month",xaxt="n",
          main="Seasonal deviation plot: Demand for 1/4 CLEAR 96X130")

ggseasonplot(A, year.labels=TRUE, continuous=TRUE)

ets(A)

## ETS(M,A,M) 
## 
## Call:
##  ets(y = A) 
## 
##   Smoothing parameters:
##     alpha = 0.001 
##     beta  = 8e-04 
##     gamma = 1e-04 
## 
##   Initial states:
##     l = 1243.9225 
##     b = 24.0393 
##     s=1.1264 1.0676 1.215 1.0808 1.1453 1.1002
##            0.9214 0.9136 0.9221 0.8216 0.8366 0.8494
## 
##   sigma:  0.1157
## 
##      AIC     AICc      BIC 
## 893.6213 908.9213 928.6489

fit1 <- ets(A)
fit2 <- ets(A,model= "MAM")
deviance <- 2*c(logLik(fit1) - logLik(fit2))
df <- attributes(logLik(fit1))$df - attributes(logLik(fit2))$df 
#P value
1-pchisq(deviance,df)

## [1] 1

ses(A)

##          Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## Nov 2016       2506.998 2077.698 2936.298 1850.440 3163.556
## Dec 2016       2506.998 1990.035 3023.961 1716.371 3297.625
## Jan 2017       2506.998 1915.218 3098.778 1601.949 3412.047
## Feb 2017       2506.998 1848.852 3165.144 1500.451 3513.545
## Mar 2017       2506.998 1788.591 3225.405 1408.290 3605.706
## Apr 2017       2506.998 1733.008 3280.988 1323.282 3690.714
## May 2017       2506.998 1681.157 3332.839 1243.984 3770.012
## Jun 2017       2506.998 1632.375 3381.621 1169.378 3844.618
## Jul 2017       2506.998 1586.174 3427.822 1098.719 3915.277
## Aug 2017       2506.998 1542.182 3471.814 1031.439 3982.557

holt(A)

##          Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## Nov 2016       2553.267 2126.425 2980.110 1900.468 3206.067
## Dec 2016       2577.451 2068.754 3086.149 1799.465 3355.437
## Jan 2017       2601.635 2022.518 3180.751 1715.953 3487.316
## Feb 2017       2625.818 1983.945 3267.692 1644.157 3607.479
## Mar 2017       2650.002 1950.966 3349.038 1580.919 3719.085
## Apr 2017       2674.185 1922.305 3426.066 1524.283 3824.088
## May 2017       2698.369 1897.107 3499.631 1472.944 3923.794
## Jun 2017       2722.553 1874.766 3570.339 1425.975 4019.130
## Jul 2017       2746.736 1854.836 3638.636 1382.692 4110.780
## Aug 2017       2770.920 1836.974 3704.865 1342.574 4199.266

train = A[1:48]
test = A[49:58]
arma_fit <- auto.arima(train)
arma_forecast <- forecast(arma_fit, h = 5)
arma_fit_accuracy <- accuracy(arma_forecast, test)
arma_fit; arma_forecast; arma_fit_accuracy

## Series: train 
## ARIMA(0,1,1)                    
## 
## Coefficients:
##           ma1
##       -0.3077
## s.e.   0.1423
## 
## sigma^2 estimated as 109140:  log likelihood=-338.84
## AIC=681.69   AICc=681.96   BIC=685.39

##    Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## 49       2709.067 2285.690 3132.444 2061.568 3356.566
## 50       2709.067 2194.132 3224.003 1921.541 3496.593
## 51       2709.067 2116.556 3301.578 1802.900 3615.235
## 52       2709.067 2048.022 3370.112 1698.086 3720.048
## 53       2709.067 1985.955 3432.179 1603.163 3814.972

##                      ME     RMSE      MAE         MPE     MAPE      MASE
## Training set   48.71659 323.4069 257.9368   0.6031599 14.31161 0.9005372
## Test set     -546.26716 589.5893 546.2672 -26.6337271 26.63373        NA
##                     ACF1
## Training set -0.02110112
## Test set              NA

plot(arma_forecast, ylim=c(0,3500))
lines(A[1:58])

Afit1 <- meanf(A,h=3)
Afit2 <- rwf(A,h=3)
Afit3 <- snaive(A,h=3)

plot(Afit1, plot.conf=FALSE,
  main="Forecasts for 1/4 CLEAR 96X130")
lines(Afit2$mean,col=2)
lines(Afit3$mean,col=3)
lines(A)
legend("topleft", lty=1, col=c(4,2,3),
  legend=c("Mean method","Naive method","Seasonal naive method"))

cycle(A)

##      Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
## 2012   1   2   3   4   5   6   7   8   9  10  11  12
## 2013   1   2   3   4   5   6   7   8   9  10  11  12
## 2014   1   2   3   4   5   6   7   8   9  10  11  12
## 2015   1   2   3   4   5   6   7   8   9  10  11  12
## 2016   1   2   3   4   5   6   7   8   9  10

boxplot(A~cycle(A))

3/32 CLEAR 72 X 84

B <- ts(c(468,538,687,932,1162,1067,1131,1201,1183,1082,793,505,498,382,401,752,955,878,999,1053,1007,1283,753,396,397,381,422,742,1028,1168,1216,1179,1329,1161,575,622,347,296,566,820,919,1379,1159,1132,1119,1238,904,699,512,471,744,789,943,1341,973,1255,1146,1045
), 
        start = c(2012, 1), frequency = 12)

plot(B, ylab="Demand for 3/32 CLEAR 72 X 84")

seasonplot(B,ylab="Demand for 3/32 CLEAR 72 X 84", xlab="Month",
           main="Seasonal plot: SS02CLR7284 Demand",
           year.labels=TRUE, year.labels.left=TRUE, col=1:20, pch=19)

what does the $df sign mean

monthplot(B,ylab="Demand for 3/32 CLEAR 72 X 84",xlab="Month",xaxt="n",
          main="Seasonal deviation plot: Demand for 3/32 CLEAR 72 X 84")

ggseasonplot(B, col =rainbow(12), year.labels=TRUE )

ets(B)

## ETS(M,N,A) 
## 
## Call:
##  ets(y = B) 
## 
##   Smoothing parameters:
##     alpha = 0.1999 
##     gamma = 0.0016 
## 
##   Initial states:
##     l = 923.9891 
##     s=-292.0386 -96.8391 334.8276 331.2985 286.2546 268.5247
##            302.8817 126.1503 -78.3819 -299.3737 -452.0562 -431.2479
## 
##   sigma:  0.1271
## 
##      AIC     AICc      BIC 
## 802.9122 814.3407 833.8188

fitB1 <- ets(B)
fitB2 <- ets(B,model= "MNA")
Bdeviance <- 2*c(logLik(fitB1) - logLik(fitB2))
Bdf <- attributes(logLik(fitB1))$df - attributes(logLik(fitB2))$df 
#P value
1-pchisq(Bdeviance,Bdf)

## [1] 1

ses(B)

##          Point Forecast    Lo 80    Hi 80      Lo 95    Hi 95
## Nov 2016        1045.01 762.6428 1327.377  613.16665 1476.854
## Dec 2016        1045.01 645.7025 1444.318  434.32177 1655.698
## Jan 2017        1045.01 555.9683 1534.052  297.08517 1792.935
## Feb 2017        1045.01 480.3179 1609.702  181.38798 1908.632
## Mar 2017        1045.01 413.6682 1676.352   79.45604 2010.564
## Apr 2017        1045.01 353.4120 1736.608  -12.69785 2102.718
## May 2017        1045.01 298.0006 1792.020  -97.44234 2187.463
## Jun 2017        1045.01 246.4248 1843.595 -176.32075 2266.341
## Jul 2017        1045.01 197.9836 1892.037 -250.40509 2340.425
## Aug 2017        1045.01 152.1668 1937.853 -320.47588 2410.496

exp <- ses(B[1:48], 5, initial="simple")
exp_accuracy = accuracy(exp, B[49:58])
exp; exp_accuracy

##    Point Forecast     Lo 80     Hi 80      Lo 95    Hi 95
## 49            699 419.49276  978.5072  271.53057 1126.469
## 50            699 303.71708 1094.2829   94.46694 1303.533
## 51            699 214.87927 1183.1207  -41.39877 1439.399
## 52            699 139.98553 1258.0145 -155.93886 1553.939
## 53            699  74.00282 1323.9972 -256.85070 1654.851

##                   ME     RMSE      MAE       MPE     MAPE      MASE
## Training set  4.8125 218.1007 165.4792 -3.668359 22.49350 0.9791667
## Test set     -7.2000 176.9825 158.8000 -8.320196 25.65224        NA
##                   ACF1
## Training set 0.2305678
## Test set            NA

plot(exp, ylim=c(0,1400))
lines(B[1:1400])

holt(B)

##          Point Forecast    Lo 80    Hi 80      Lo 95    Hi 95
## Nov 2016       1063.694 780.1760 1347.212  630.09064 1497.298
## Dec 2016       1078.575 683.9211 1473.230  475.00380 1682.147
## Jan 2017       1093.457 612.6895 1574.224  358.18679 1828.726
## Feb 2017       1108.338 554.6795 1661.996  261.59044 1955.085
## Mar 2017       1123.219 505.1950 1741.243  178.03283 2068.405
## Apr 2017       1138.100 461.7975 1814.403  103.78444 2172.416
## May 2017       1152.982 423.0276 1882.936   36.61331 2269.350
## Jun 2017       1167.863 387.9297 1947.796  -24.94197 2360.668
## Jul 2017       1182.744 355.8377 2009.651  -81.90010 2447.388
## Aug 2017       1197.625 326.2654 2068.985 -135.00476 2530.256

Btrain = B[1:48]
Btest = B[49:58]
arma_fitB <- auto.arima(Btrain)
arma_forecastB <- forecast(arma_fitB, h = 5)
arma_fit_accuracyB <- accuracy(arma_forecastB, Btest)
arma_fitB; arma_forecastB; arma_fit_accuracyB

## Series: Btrain 
## ARIMA(4,0,1) with non-zero mean 
## 
## Coefficients:
##          ar1      ar2     ar3      ar4      ma1  intercept
##       1.2689  -0.5530  0.1912  -0.3419  -0.7196   836.2356
## s.e.  0.1674   0.2375  0.2280   0.1514   0.1400    15.1846
## 
## sigma^2 estimated as 23791:  log likelihood=-308.52
## AIC=631.04   AICc=633.84   BIC=644.14

##    Point Forecast    Lo 80    Hi 80     Lo 95     Hi 95
## 49       563.0191 365.3486 760.6896 260.70821  865.3300
## 50       441.0354 215.5021 666.5686  96.11204  785.9587
## 51       436.4414 209.1167 663.7661  88.77823  784.1045
## 52       542.1569 314.4088 769.9051 193.84619  890.4677
## 53       702.0119 469.8689 934.1548 346.97989 1057.0438

##                       ME     RMSE      MAE      MPE     MAPE      MASE
## Training set   0.6824324 144.2812 112.3856 -4.06676 15.36994 0.6650034
## Test set     154.8670725 208.3746 175.2747 18.91536 22.90123        NA
##                     ACF1
## Training set -0.01456141
## Test set              NA

plot(arma_forecastB, ylim=c(0,2000))
lines(B[1:58])

B2 <- window(B, start = 2012.1, end = 2016-.10)
Bfit1 <- meanf(B2,h=3)
Bfit2 <- rwf(B2,h=3)
Bfit3 <- snaive(B2,h=3)

plot(Bfit1, plot.conf=FALSE,
  main="Forecasts for 3/32 CLEAR 72 X 84")
lines(Bfit2$mean,col=2)
lines(Bfit3$mean,col=3)
lines(B)
legend("topleft", lty=1, col=c(4,2,3),
  legend=c("Mean method","Naive method","Seasonal naive method"))

B3 <- window(B, start = 2016-.10)
accuracy(Bfit1, B3)

##                         ME     RMSE      MAE       MPE     MAPE     MASE
## Training set  4.294196e-14 315.5095 276.2528 -20.40474 44.02496 1.918422
## Test set     -3.104222e+02 325.8990 310.4222 -59.89947 59.89947 2.155710
##                     ACF1 Theil's U
## Training set  0.74896028        NA
## Test set     -0.08016487  3.348802

accuracy(Bfit2, B3)

##                       ME     RMSE      MAE        MPE     MAPE     MASE
## Training set    4.931818 224.3239 170.8864  -4.123938 23.08320 1.186711
## Test set     -343.333333 357.3877 343.3333 -65.940723 65.94072 2.384259
##                     ACF1 Theil's U
## Training set  0.18588381        NA
## Test set     -0.08016487  3.633793

accuracy(Bfit3, B3)

##                     ME     RMSE MAE       MPE     MAPE      MASE
## Training set  -3.69697 168.6032 144 -4.025599 18.23994 1.0000000
## Test set     139.00000 145.8069 139 26.799096 26.79910 0.9652778
##                    ACF1 Theil's U
## Training set  0.2487182        NA
## Test set     -0.1162311   1.48749

cycle(B)

##      Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
## 2012   1   2   3   4   5   6   7   8   9  10  11  12
## 2013   1   2   3   4   5   6   7   8   9  10  11  12
## 2014   1   2   3   4   5   6   7   8   9  10  11  12
## 2015   1   2   3   4   5   6   7   8   9  10  11  12
## 2016   1   2   3   4   5   6   7   8   9  10

boxplot(B~cycle(B))

fit.B <- tslm(B~cycle(B))
plot(fit.B)

BBB <- forecast(fit.B,  h=58,level=c(80,95))
plot(BBB, ylab = "Demand", xlab = "Time")
lines(fitted(fit.B), col = "blue")

summary(fit.B)

## 
## Call:
## tslm(formula = B ~ cycle(B))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -714.89 -194.86   42.55  211.22  529.06 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   589.00      77.47   7.603 3.48e-10 ***
## cycle(B)       43.49      10.80   4.028 0.000171 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 278 on 56 degrees of freedom
## Multiple R-squared:  0.2246, Adjusted R-squared:  0.2108 
## F-statistic: 16.22 on 1 and 56 DF,  p-value: 0.0001714

1/4 SOLARBAN 60 TC 100 X 144

C <- ts(c(104,105,99,94,116,113,140,134,97,95,115,61,185,122,125,126,36,85,249,198,185,284,328,392,314,235,120,206,198,332,200,458,437,441,338,330,316,300,294,321,238,415,652,424,469,414,637,612,605,487,478,575,682,784,653,903,612,590
), 
        start = c(2012, 1), frequency = 12)
plot(C)

seasonplot(C,ylab="Demand", xlab="Month",
           main="Seasonal plot: 1/4 SOLARBAN 60 TC 100 X 144 Demand",
           year.labels=TRUE, year.labels.left=TRUE, col=1:20, pch=19)

monthplot(C,ylab="Demand",xlab="Month",xaxt="n",
          main="1/4 SOLARBAN 60 TC 100 X 144")

ets(C)

## ETS(M,A,N) 
## 
## Call:
##  ets(y = C) 
## 
##   Smoothing parameters:
##     alpha = 0.2046 
##     beta  = 1e-04 
## 
##   Initial states:
##     l = 111.6516 
##     b = 8.5027 
## 
##   sigma:  0.3268
## 
##      AIC     AICc      BIC 
## 766.5066 767.6605 776.8089

Cfit1 <- ets(C)
Cfit2 <- ets(C,model= "MAN")
Cdeviance <- 2*c(logLik(Cfit1) - logLik(Cfit2))
Cdf <- attributes(logLik(Cfit1))$Cdf - attributes(logLik(Cfit2))$Cdf 
#P value
1-pchisq(Cdeviance,Cdf)

## numeric(0)

ses(C)

##          Point Forecast    Lo 80    Hi 80    Lo 95     Hi 95
## Nov 2016       640.3345 515.5047 765.1644 449.4237  831.2454
## Dec 2016       640.3345 498.6880 781.9811 423.7049  856.9642
## Jan 2017       640.3345 483.6662 797.0029 400.7310  879.9381
## Feb 2017       640.3345 469.9638 810.7053 379.7749  900.8942
## Mar 2017       640.3345 457.2842 823.3849 360.3832  920.2859
## Apr 2017       640.3345 445.4277 835.2414 342.2503  938.4188
## May 2017       640.3345 434.2523 846.4168 325.1589  955.5102
## Jun 2017       640.3345 423.6524 857.0167 308.9479  971.7212
## Jul 2017       640.3345 413.5475 867.1216 293.4937  987.1754
## Aug 2017       640.3345 403.8740 876.7951 278.6993 1001.9698

holt(C)

##          Point Forecast    Lo 80    Hi 80    Lo 95     Hi 95
## Nov 2016       695.7371 574.2254 817.2487 509.9010  881.5731
## Dec 2016       706.2466 580.2455 832.2477 513.5446  898.9487
## Jan 2017       716.7562 586.4171 847.0953 517.4198  916.0927
## Feb 2017       727.2658 592.7255 861.8061 521.5042  933.0274
## Mar 2017       737.7754 599.1582 876.3925 525.7788  949.7720
## Apr 2017       748.2850 605.7047 890.8653 530.2272  966.3427
## May 2017       758.7945 612.3555 905.2336 534.8354  982.7537
## Jun 2017       769.3041 619.1028 919.5055 539.5910  999.0173
## Jul 2017       779.8137 625.9394 933.6881 544.4832 1015.1442
## Aug 2017       790.3233 632.8590 947.7876 549.5025 1031.1441

Ctrain = C[1:48]
Ctest = C[49:58]
Carma_fit <- auto.arima(Ctrain)
Carma_forecast <- forecast(Carma_fit, h = 5)
Carma_fit_accuracy <- accuracy(Carma_forecast, Ctest)
Carma_fit; Carma_forecast; Carma_fit_accuracy

## Series: Ctrain 
## ARIMA(0,1,1)                    
## 
## Coefficients:
##           ma1
##       -0.3906
## s.e.   0.1628
## 
## sigma^2 estimated as 8295:  log likelihood=-278.32
## AIC=560.63   AICc=560.91   BIC=564.34

##    Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## 49        590.993 474.2735 707.7125 412.4859 769.5000
## 50        590.993 454.3110 727.6750 381.9559 800.0300
## 51        590.993 436.9135 745.0725 355.3487 826.6372
## 52        590.993 421.2902 760.6957 331.4550 850.5309
## 53        590.993 406.9888 774.9972 309.5828 872.4032

##                     ME     RMSE     MAE       MPE     MAPE      MASE
## Training set  16.65187 89.15907 60.7501 -4.912544 27.58987 0.9180883
## Test set     -25.59298 80.39428 67.5986 -6.422906 12.68665        NA
##                      ACF1
## Training set -0.003465877
## Test set               NA

plot(Carma_forecast, ylim=c(0,1000))
lines(C[1:58])

Cfit1 <- meanf(C,h=3)
Cfit2 <- rwf(C,h=3)
Cfit3 <- snaive(C,h=3)

plot(Cfit1, plot.conf=FALSE,
  main="Forecasts for 1/4 SOLARBAN 60 TC 100 X 144")
lines(Cfit2$mean,col=2)
lines(Cfit3$mean,col=3)
lines(C)
legend("topleft", lty=1, col=c(4,2,3),
  legend=c("Mean method","Naive method","Seasonal naive method"))

cycle(C)

##      Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
## 2012   1   2   3   4   5   6   7   8   9  10  11  12
## 2013   1   2   3   4   5   6   7   8   9  10  11  12
## 2014   1   2   3   4   5   6   7   8   9  10  11  12
## 2015   1   2   3   4   5   6   7   8   9  10  11  12
## 2016   1   2   3   4   5   6   7   8   9  10

boxplot(C~cycle(C))

3/16 CLEAR 96X130

D <- ts(c(278,366,381,393,498,441,381,534,453,534,391,288,302,319,384,372,381,371,407,365,532,445,367,266,250,266,168,323,360,329,432,395,471,549,295,302,219,179,321,389,299,488,449,509,499,561,466,373,262,337,408,437,420,704,492,726,768,519
), 
        start = c(2012, 1), frequency = 12)
plot(D)

seasonplot(D,ylab="Demand", xlab="Month",
           main="Seasonal plot: 3/16 CLEAR 96X130 Demand",
           year.labels=TRUE, year.labels.left=TRUE, col=1:20, pch=19)

monthplot(C,ylab="Demand",xlab="Month",xaxt="n",
          main="1/4 SOLARBAN 60 TC 100 X 144")

ets(D)

## ETS(M,A,A) 
## 
## Call:
##  ets(y = D) 
## 
##   Smoothing parameters:
##     alpha = 0.0377 
##     beta  = 0.0377 
##     gamma = 1e-04 
## 
##   Initial states:
##     l = 448.7472 
##     b = -3.5758 
##     s=-70.4183 2.9142 145.0808 112.7078 74.4152 40.423
##            35.6681 -16.7346 -13.6575 -69.228 -109.8366 -131.3342
## 
##   sigma:  0.1496
## 
##      AIC     AICc      BIC 
## 738.4683 753.7683 773.4958

Dfit1 <- ets(D)
Dfit2 <- ets(D,model= "MAA")
Ddeviance <- 2*c(logLik(Dfit1) - logLik(Dfit2))
Ddf <- attributes(logLik(Dfit1))$Ddf - attributes(logLik(Dfit2))$Ddf 
#P value
1-pchisq(Ddeviance,Ddf)

## numeric(0)

ses(D)

##          Point Forecast    Lo 80    Hi 80    Lo 95     Hi 95
## Nov 2016       591.5978 463.2344 719.9611 395.2829  787.9126
## Dec 2016       591.5978 438.0055 745.1901 356.6986  826.4969
## Jan 2017       591.5978 416.3721 766.8234 323.6132  859.5823
## Feb 2017       591.5978 397.1306 786.0649 294.1859  889.0096
## Mar 2017       591.5978 379.6286 803.5669 267.4190  915.7765
## Apr 2017       591.5978 363.4655 819.7301 242.6995  940.4960
## May 2017       591.5978 348.3740 834.8215 219.6192  963.5763
## Jun 2017       591.5978 334.1658 849.0297 197.8895  985.3060
## Jul 2017       591.5978 320.7017 862.4938 177.2980 1005.8975
## Aug 2017       591.5978 307.8759 875.3196 157.6826 1025.5129

holt(D)

##          Point Forecast    Lo 80    Hi 80    Lo 95     Hi 95
## Nov 2016       601.1418 473.0810 729.2026 405.2897  796.9939
## Dec 2016       605.9664 453.6793 758.2534 373.0634  838.8693
## Jan 2017       610.7909 437.6285 783.9534 345.9618  875.6201
## Feb 2017       615.6155 423.8312 807.3999 322.3066  908.9244
## Mar 2017       620.4401 411.6834 829.1968 301.1743  939.7059
## Apr 2017       625.2647 400.8107 849.7186 281.9919  968.5374
## May 2017       630.0892 390.9618 869.2166 264.3754  995.8030
## Jun 2017       634.9138 381.9585 887.8691 248.0521 1021.7755
## Jul 2017       639.7384 373.6690 905.8078 232.8203 1046.6564
## Aug 2017       644.5629 365.9923 923.1336 218.5259 1070.6000

Dtrain = D[1:48]
Dtest = D[49:58]
Darma_fit <- auto.arima(Dtrain)
Darma_forecast <- forecast(Darma_fit, h = 5)
Darma_fit_accuracy <- accuracy(Darma_forecast, Dtest)
Darma_fit; Darma_forecast; Darma_fit_accuracy

## Series: Dtrain 
## ARIMA(1,0,0) with non-zero mean 
## 
## Coefficients:
##          ar1  intercept
##       0.5550   379.3081
## s.e.  0.1188    25.1007
## 
## sigma^2 estimated as 6557:  log likelihood=-278.19
## AIC=562.38   AICc=562.93   BIC=568

##    Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## 49       375.8072 272.0311 479.5833 217.0954 534.5190
## 50       377.3652 258.6785 496.0518 195.8496 558.8807
## 51       378.2298 255.3145 501.1452 190.2470 566.2126
## 52       378.7097 254.5209 502.8985 188.7793 568.6401
## 53       378.9760 254.3976 503.5544 188.4497 569.5022

##                     ME     RMSE      MAE       MPE     MAPE      MASE
## Training set  1.526237 79.27194 62.93519 -4.813821 18.48928 0.8661652
## Test set     -5.017570 64.10683 56.65139 -5.002536 17.16373        NA
##                     ACF1
## Training set -0.03122273
## Test set              NA

plot(Darma_forecast, ylim=c(0,1000))
lines(D[1:58])

Dfit1 <- meanf(D,h=3)
Dfit2 <- rwf(D,h=3)
Dfit3 <- snaive(D,h=3)

plot(Dfit1, plot.conf=FALSE,
  main="Forecasts for 3/16 CLEAR 96X130")
lines(Dfit2$mean,col=2)
lines(Dfit3$mean,col=3)
lines(D)
legend("topleft", lty=1, col=c(4,2,3),
  legend=c("Mean method","Naive method","Seasonal naive method"))

cycle(D)

##      Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
## 2012   1   2   3   4   5   6   7   8   9  10  11  12
## 2013   1   2   3   4   5   6   7   8   9  10  11  12
## 2014   1   2   3   4   5   6   7   8   9  10  11  12
## 2015   1   2   3   4   5   6   7   8   9  10  11  12
## 2016   1   2   3   4   5   6   7   8   9  10

boxplot(D~cycle(D))

3/8 CLEAR 96 X 130

E <- ts(c(196,238,282,250,309,269,215,314,185,242,188,151,196,292,244,199,320,384,272,263,225,215,265,219,202,236,274,234,275,298,320,362,426,408,244,280,211,253,410,449,408,450,441,462,470,396,468,544,293,415,467,402,490,598,501,380,529,439
), start = c(2012, 1), frequency = 12)


plot(E)

seasonplot(E,ylab="Demand for 3/8 CLEAR 96 X 130
", xlab="Month",
           main="Seasonal plot: 3/8 CLEAR 96 X 130
 Demand",
           year.labels=TRUE, year.labels.left=TRUE, col=1:20, pch=19)

monthplot(E,ylab="Demand for 3/8 CLEAR 96 X 130",xlab="Month",xaxt="n",
          main="Seasonal deviation plot: Demand for 3/8 CLEAR 96 X 130")

ets(E)

## ETS(M,N,N) 
## 
## Call:
##  ets(y = E) 
## 
##   Smoothing parameters:
##     alpha = 0.4002 
## 
##   Initial states:
##     l = 225.0354 
## 
##   sigma:  0.2254
## 
##      AIC     AICc      BIC 
## 731.9426 732.3870 738.1239

fitE1 <- ets(E)
fitE2 <- ets(E,model= "MNN")
Edeviance <- 2*c(logLik(fitE1) - logLik(fitE2))
Edf <- attributes(logLik(fitE1))$df - attributes(logLik(fitE2))$df 
#P value
1-pchisq(Edeviance,Edf)

## [1] 1

ses(E)

##          Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## Nov 2016       466.9117 376.8411 556.9824 329.1605 604.6629
## Dec 2016       466.9117 370.5904 563.2330 319.6010 614.2224
## Jan 2017       466.9117 364.7214 569.1020 310.6252 623.1983
## Feb 2017       466.9117 359.1717 574.6518 302.1375 631.6859
## Mar 2017       466.9117 353.8941 579.9293 294.0662 639.7573
## Apr 2017       466.9117 348.8522 584.9712 286.3553 647.4682
## May 2017       466.9117 344.0170 589.8064 278.9605 654.8630
## Jun 2017       466.9117 339.3650 594.4585 271.8458 661.9777
## Jul 2017       466.9117 334.8767 598.9467 264.9816 668.8418
## Aug 2017       466.9117 330.5361 603.2873 258.3432 675.4802

holt(E)

##          Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## Nov 2016       485.1267 396.6159 573.6375 349.7611 620.4923
## Dec 2016       489.4753 398.7988 580.1518 350.7975 628.1531
## Jan 2017       493.8239 401.0303 586.6175 351.9083 635.7395
## Feb 2017       498.1725 403.3072 593.0378 353.0886 643.2564
## Mar 2017       502.5211 405.6266 599.4155 354.3338 650.7084
## Apr 2017       506.8697 407.9859 605.7535 355.6399 658.0994
## May 2017       511.2183 410.3826 612.0539 357.0035 665.4330
## Jun 2017       515.5668 412.8148 618.3189 358.4211 672.7126
## Jul 2017       519.9154 415.2803 624.5506 359.8898 679.9410
## Aug 2017       524.2640 417.7775 630.7505 361.4070 687.1211

Dtrain = E[1:48]
Dtest = B[49:58]
arma_fitE <- auto.arima(Dtrain)
arma_forecastE <- forecast(arma_fitE, h = 5)
arma_fit_accuracyE <- accuracy(arma_forecastE, Dtest)
arma_fitE; arma_forecastE; arma_fit_accuracyE

## Series: Dtrain 
## ARIMA(0,1,0)                    
## 
## sigma^2 estimated as 4139:  log likelihood=-262.4
## AIC=526.8   AICc=526.89   BIC=528.65

##    Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## 49            544 461.5562 626.4438 417.9130 670.0870
## 50            544 427.4068 660.5932 365.6861 722.3139
## 51            544 401.2031 686.7969 325.6110 762.3890
## 52            544 379.1123 708.8877 291.8261 796.1739
## 53            544 359.6500 728.3500 262.0610 825.9390

##                      ME      RMSE       MAE        MPE     MAPE      MASE
## Training set   7.254083  63.65762  52.50408 -0.4853897 18.92227 0.9792428
## Test set     147.800000 230.46865 189.80000 15.6993035 24.39888        NA
##                    ACF1
## Training set -0.1744284
## Test set             NA

plot(arma_forecastE, ylim = c(0,600))
lines(E[1:58])

E2 <- window(E, start = 2012.1, end = 2016-.11)
Efit1 <- meanf(E2,h=3)
Efit2 <- rwf(E2,h=3)
Efit3 <- snaive(E2,h=3)

plot(Efit1, plot.conf=FALSE,
  main="Forecasts for 3/8 CLEAR 96 X 130")
lines(Efit2$mean,col=2)
lines(Efit3$mean,col=3)
lines(E)
legend("topleft", lty=1, col=c(4,2,3),
  legend=c("Mean method","Naive method","Seasonal naive method"))

E3 <- window(E, start = 2016-.10)
accuracy(Efit1, E3)

##                        ME      RMSE       MAE       MPE     MAPE      MASE
## Training set 2.652230e-14  89.79126  75.82519 -9.019279 26.81491 0.9702331
## Test set     1.178667e+02 156.19040 122.17778 23.527761 24.99913 1.5633450
##                    ACF1 Theil's U
## Training set  0.7054355        NA
## Test set     -0.4906199 0.6347406

boxplot(E~cycle(E))

FF <- ts(c(168,81,136,168,148,103,101,194,256,232,160,130,190,231,215,297,222,336,326,436,221,185,255,217,160,153,198,230,263,307,319,330,433,384,275,358,248,307,364,390,384,429,622,498,485,423,291,320,325,465,440,416,475,509,642,677,437,414
), start = c(2012,1), frequency = 12)
plot(FF)

seasonplot(FF, ylab = "Demand for 1/2 CLEAR 102X130", xlab="Month", main = "Seasonal Plot:1/2 CLEAR 102X130 Demand", year.labels = TRUE, year.labels.left = TRUE, col = 1:20, pch= 19)

monthplot(FF,ylab="Demand for 1/2 CLEAR 102X130",xlab="Month",xaxt="n",
          main="Seasonal deviation plot: Demand for 1/2 CLEAR 102X130")

ggseasonplot(FF,year.labels=TRUE, continuous=TRUE)

ets(FF)

## ETS(M,A,N) 
## 
## Call:
##  ets(y = FF) 
## 
##   Smoothing parameters:
##     alpha = 0.6804 
##     beta  = 1e-04 
## 
##   Initial states:
##     l = 109.641 
##     b = 11.8682 
## 
##   sigma:  0.2475
## 
##      AIC     AICc      BIC 
## 741.7355 742.8893 752.0377

fitF1 <- ets(FF)
fitF2 <- ets(FF,model= "MAN")
Fdeviance <- 2*c(logLik(fitF1) - logLik(fitF2))
Fdf <- attributes(logLik(fitF1))$Fdf - attributes(logLik(fitF2))$Fdf 
#P value
1-pchisq(Fdeviance,Fdf)

## numeric(0)

ses(FF)

##          Point Forecast    Lo 80    Hi 80     Lo 95    Hi 95
## Nov 2016       428.6912 327.6339 529.7485 274.13734 583.2450
## Dec 2016       428.6912 299.8365 557.5458 231.62497 625.7574
## Jan 2017       428.6912 277.0519 580.3304 196.77892 660.6034
## Feb 2017       428.6912 257.2694 600.1129 166.52425 690.8581
## Mar 2017       428.6912 239.5448 617.8376 139.41674 717.9656
## Apr 2017       428.6912 223.3444 634.0380 114.64038 742.7420
## May 2017       428.6912 208.3318 649.0505  91.68064 765.7017
## Jun 2017       428.6912 194.2788 663.1036  70.18831 787.1940
## Jul 2017       428.6912 181.0218 676.3606  49.91353 807.4688
## Aug 2017       428.6912 168.4392 688.9431  30.67019 826.7122

holt(FF)

##          Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## Nov 2016       439.3646 338.4636 540.2656 285.0498 593.6794
## Dec 2016       444.5563 317.8775 571.2351 250.8178 638.2948
## Jan 2017       449.7480 301.7090 597.7871 223.3418 676.1542
## Feb 2017       454.9397 288.2509 621.6285 200.0112 709.8682
## Mar 2017       460.1314 276.6749 643.5879 179.5589 740.7040
## Apr 2017       465.3231 266.5042 664.1421 161.2558 769.3905
## May 2017       470.5148 257.4345 683.5951 144.6366 796.3930
## Jun 2017       475.7065 249.2579 702.1552 129.3832 822.0299
## Jul 2017       480.8982 241.8243 719.9721 115.2663 846.5302
## Aug 2017       486.0899 235.0218 737.1580 102.1144 870.0655

Ftrain = FF[1:48]
Ftest = FF[49:58]
arma_fitF <- auto.arima(Ftrain)
arma_forecastF <- forecast(arma_fitF, h = 5)
arma_fit_accuracyF <- accuracy(arma_forecastF, Ftest)
arma_fitF; arma_forecastF; arma_fit_accuracyF

## Series: Ftrain 
## ARIMA(0,1,0)                    
## 
## sigma^2 estimated as 5703:  log likelihood=-269.94
## AIC=541.87   AICc=541.96   BIC=543.72

##    Point Forecast    Lo 80    Hi 80     Lo 95    Hi 95
## 49            320 223.2208 416.7792 171.98895 468.0110
## 50            320 183.1335 456.8665 110.68077 529.3192
## 51            320 152.3735 487.6265  63.63735 576.3627
## 52            320 126.4416 513.5584  23.97791 616.0221
## 53            320 103.5951 536.4049 -10.96276 650.9628

##                      ME      RMSE       MAE       MPE     MAPE      MASE
## Training set   3.170167  74.72645  58.96183 -3.013949 24.05579 0.9792248
## Test set     104.200000 117.21007 104.20000 23.140497 23.14050        NA
##                    ACF1
## Training set -0.2051327
## Test set             NA

plot(arma_forecastF, ylim=c(0,800))
lines(FF[1:58])

F2 <- window(FF, start = 2012.1, end = 2016.11)
Ffit1 <- meanf(F2,h=3)
Ffit2 <- rwf(F2,h=3)
Ffit3 <- snaive(F2,h=3)
plot(Ffit1, plot.conf=FALSE,
  main="Forecasts for 1/2 CLEAR 102X130")
lines(Ffit2$mean,col=2)
lines(Ffit3$mean,col=3)
lines(FF)
legend("topleft", lty=1, col=c(4,2,3),
  legend=c("Mean method","Naive method","Seasonal naive method"))

F3 <- window(FF, start = 2016-.11)
accuracy(Ffit1, F3)

##                        ME     RMSE       MAE       MPE     MAPE      MASE
## Training set 1.892699e-14 114.4433  93.95139 -19.10707 40.90001 0.8773671
## Test set     1.578333e+02 159.6817 157.83333  35.38417 35.38417 1.4739300
##                    ACF1 Theil's U
## Training set  0.7362556        NA
## Test set     -0.4347469  3.570083

accuracy(Ffit2, F3)

##                     ME     RMSE     MAE        MPE      MAPE      MASE
## Training set   7.00000 76.77890 60.2766 -0.9820706 22.093091 0.5628943
## Test set     -21.33333 32.28002 28.0000 -5.1184671  6.521976 0.2614786
##                    ACF1 Theil's U
## Training set -0.1883339        NA
## Test set     -0.4347469 0.7497556

accuracy(Ffit3, F3)

##                    ME      RMSE       MAE      MPE     MAPE      MASE
## Training set 81.86111 130.51724 107.08333 20.58936 31.94081 1.0000000
## Test set     64.33333  70.07853  64.33333 14.22687 14.22687 0.6007782
##                    ACF1 Theil's U
## Training set  0.4249932        NA
## Test set     -0.6342926  1.491177

cycle(FF)

##      Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
## 2012   1   2   3   4   5   6   7   8   9  10  11  12
## 2013   1   2   3   4   5   6   7   8   9  10  11  12
## 2014   1   2   3   4   5   6   7   8   9  10  11  12
## 2015   1   2   3   4   5   6   7   8   9  10  11  12
## 2016   1   2   3   4   5   6   7   8   9  10

boxplot(FF~cycle(FF))

add 15 months

1/8 ENEG ADV LOW-E 72 X 84

G <- ts(c(243,209,163,181,123,121,268,275,360,149,174,346,395,234,243,
 185,218,183,281,331,434,428,375,473,536,333,292,259,193,281,299,190,292,365,522,457,398), start = c(2013,10), frequency = 12)
plot(G)

seasonplot(G, ylab = "Demand for 1/8 ENEG ADV LOW-E 72 X 84", xlab="Month", main = "Seasonal Plot:1/8 ENEG ADV LOW-E 72 X 84 Demand", year.labels = TRUE, year.labels.left = TRUE, col = 1:20, pch= 19)

monthplot(G,ylab="Demand for 1/8 ENEG ADV LOW-E 72 X 84",xlab="Month",xaxt="n",
          main="Seasonal deviation plot: Demand for 1/8 ENEG ADV LOW-E 72 X 84")

ets(G)

## ETS(A,N,N) 
## 
## Call:
##  ets(y = G) 
## 
##   Smoothing parameters:
##     alpha = 0.8439 
## 
##   Initial states:
##     l = 236.7237 
## 
##   sigma:  89.5984
## 
##      AIC     AICc      BIC 
## 472.2589 472.9862 477.0917

fitG1 <- ets(G)
fitG2 <- ets(G,model= "ANN")
Gdeviance <- 2*c(logLik(fitG1) - logLik(fitG2))
Gdf <- attributes(logLik(fitG1))$df - attributes(logLik(fitG2))$df 
#P value
1-pchisq(Gdeviance,Gdf)

## [1] 1

ses(G)

##          Point Forecast     Lo 80    Hi 80      Lo 95    Hi 95
## Nov 2016       408.1542 293.32921 522.9791 232.544565 583.7637
## Dec 2016       408.1542 257.91581 558.3925 178.384436 637.9239
## Jan 2017       408.1542 229.38516 586.9231 134.750572 681.5577
## Feb 2017       408.1542 204.81909 611.4892  97.180017 719.1283
## Mar 2017       408.1542 182.91664 633.3917  63.683097 752.6252
## Apr 2017       408.1542 162.96294 653.3454  33.166549 783.1418
## May 2017       408.1542 144.51515 671.7932   4.953088 811.3552
## Jun 2017       408.1542 127.27639 689.0319 -21.411318 837.7196
## Jul 2017       408.1542 111.03614 705.2722 -46.248629 862.5569
## Aug 2017       408.1542  95.63871 720.6696 -69.796978 886.1053

holt(G)

##          Point Forecast    Lo 80    Hi 80     Lo 95    Hi 95
## Nov 2016       419.1148 304.2357 533.9938 243.42239 594.8071
## Dec 2016       425.8669 278.4162 573.3175 200.36060 651.3731
## Jan 2017       432.6190 258.5836 606.6543 166.45489 698.7830
## Feb 2017       439.3711 242.2999 636.4422 137.97677 740.7653
## Mar 2017       446.1232 228.4356 663.8108 113.19871 779.0476
## Apr 2017       452.8753 216.3570 689.3935  91.15184 814.5987
## May 2017       459.6274 205.6668 713.5879  71.22814 848.0266
## Jun 2017       466.3795 196.0959 736.6630  53.01638 879.7425
## Jul 2017       473.1315 187.4524 758.8107  36.22299 910.0401
## Aug 2017       479.8836 179.5937 780.1736  20.62978 939.1375

Gtrain = G[1:32]
Gtest = G[33:37]
arma_fitG <-  auto.arima(Gtrain)
arma_forecastG <- forecast(arma_fitG, h = 5)
arma_fit_accuracyG <- accuracy(arma_forecastG, Gtest)
arma_fitG; arma_forecastG; arma_fit_accuracyG

## Series: Gtrain 
## ARIMA(0,1,0)                    
## 
## sigma^2 estimated as 8131:  log likelihood=-183.54
## AIC=369.08   AICc=369.22   BIC=370.52

##    Point Forecast     Lo 80    Hi 80      Lo 95    Hi 95
## 33            190  74.43709 305.5629   13.26179 366.7382
## 34            190  26.56937 353.4306  -59.94558 439.9456
## 35            190 -10.16082 390.1608 -116.11956 496.1196
## 36            190 -41.12581 421.1258 -163.47643 543.4764
## 37            190 -68.40651 448.4065 -205.19866 585.1987

##                      ME      RMSE       MAE       MPE     MAPE     MASE
## Training set  -1.648656  88.75406  68.10134 -6.847836 27.46662 0.968858
## Test set     216.800000 230.56713 216.80000 51.432812 51.43281 3.084351
##                     ACF1
## Training set -0.08735663
## Test set              NA

plot(arma_forecastG, ylim=c(0,800))
lines(G[1:37])

G2 <- window(G, start = 2013.10, end = 2016.11)

## Warning in window.default(x, ...): 'start' value not changed

Gfit1 <- meanf(G2,h=3)
Gfit2 <- rwf(G2,h=3)
Gfit3 <- snaive(G2,h=3)

plot(Gfit1, plot.conf=FALSE,
  main="Forecasts for 1/8 ENEG ADV LOW-E 72 X 84")
lines(Gfit2$mean,col=2)
lines(Gfit3$mean,col=3)
lines(G)
legend("topleft", lty=1, col=c(4,2,3),
  legend=c("Mean method","Naive method","Seasonal naive method"))

G3 <- window(G, start = 2016-.11)
accuracy(Gfit1, G3)

##                         ME      RMSE      MAE       MPE     MAPE      MASE
## Training set -7.849564e-15 106.46012 87.69560 -16.27989 37.17971 0.9581139
## Test set     -1.936782e+01  51.49111 37.98851 -11.94448 18.24306 0.4150415
##                    ACF1 Theil's U
## Training set  0.6347348        NA
## Test set     -0.2624394 0.8081469

accuracy(Gfit2, G3)

##                     ME     RMSE      MAE       MPE     MAPE      MASE
## Training set -1.785714 91.05022 70.14286 -7.114254 28.00452 0.7663423
## Test set     63.666667 79.55920 65.66667 21.729761 22.78239 0.7174379
##                     ACF1 Theil's U
## Training set -0.06818766        NA
## Test set     -0.26243939  1.019517

accuracy(Gfit3, G3)

##                     ME      RMSE      MAE       MPE     MAPE      MASE
## Training set 88.588235 114.87845 91.52941  25.55137 27.07530 1.0000000
## Test set     -8.333333  99.68116 85.66667 -11.10500 38.36868 0.9359469
##                     ACF1 Theil's U
## Training set  0.31672727        NA
## Test set     -0.02342665  1.285759

cycle(G)

##      Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
## 2013                                      10  11  12
## 2014   1   2   3   4   5   6   7   8   9  10  11  12
## 2015   1   2   3   4   5   6   7   8   9  10  11  12
## 2016   1   2   3   4   5   6   7   8   9  10

boxplot(G~cycle(G))

 H <- ts(c(69,140,81,50,35,56,63,180,141,122,56,153,97,81,43,91,154,63,186,89,67,121,164,259,156,185,84,122,92,96,184,144,191,201,172,224,269,151,152,149,134,212,246,165,232,382,358,336,194,232,272,288,289,239,248,284,298,398
), 
        start = c(2012, 1), frequency = 12)

plot(H)

seasonplot(H,ylab="Demand", xlab="Month",
           main="Seasonal plot: 1/4 SOLARBAN 60 100 X 144 Demand",
           year.labels=TRUE, year.labels.left=TRUE, col=1:20, pch=19)

monthplot(H,ylab="Demand",xlab="Month",xaxt="n",
          main="Seasonal deviation plot: Demand for 1/4 SOLARBAN 60 100 X 144")

ggseasonplot(H, year.labels=TRUE, continuous=TRUE)

ets(H)

## ETS(A,A,N) 
## 
## Call:
##  ets(y = H) 
## 
##   Smoothing parameters:
##     alpha = 1e-04 
##     beta  = 1e-04 
## 
##   Initial states:
##     l = 66.3181 
##     b = 3.7355 
## 
##   sigma:  55.3721
## 
##      AIC     AICc      BIC 
## 711.1385 712.2924 721.4407

Hfit1 <- ets(H)
Hfit2 <- ets(H,model= "AAN")
Hdeviance <- 2*c(logLik(Hfit1) - logLik(Hfit2))
Hdf <- attributes(logLik(Hfit1))$Hdf - attributes(logLik(Hfit2))$Hdf 
#P value
1-pchisq(Hdeviance,Hdf)

## numeric(0)

ses(H)

##          Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## Nov 2016       323.7386 248.1573 399.3200 208.1469 439.3304
## Dec 2016       323.7386 242.9927 404.4846 200.2483 447.2289
## Jan 2017       323.7386 238.1391 409.3382 192.8255 454.6518
## Feb 2017       323.7386 233.5464 413.9309 185.8015 461.6758
## Mar 2017       323.7386 229.1764 418.3009 179.1182 468.3591
## Apr 2017       323.7386 224.9997 422.4776 172.7304 474.7468
## May 2017       323.7386 220.9926 426.4847 166.6021 480.8752
## Jun 2017       323.7386 217.1360 430.3413 160.7040 486.7733
## Jul 2017       323.7386 213.4142 434.0631 155.0120 492.4653
## Aug 2017       323.7386 209.8139 437.6634 149.5058 497.9715

holt(H)

##          Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## Nov 2016       285.0992 214.1370 356.0614 176.5718 393.6265
## Dec 2016       288.8085 217.8463 359.7707 180.2811 397.3358
## Jan 2017       292.5178 221.5556 363.4800 183.9904 401.0451
## Feb 2017       296.2271 225.2649 367.1893 187.6997 404.7544
## Mar 2017       299.9364 228.9742 370.8986 191.4090 408.4638
## Apr 2017       303.6457 232.6834 374.6079 195.1183 412.1731
## May 2017       307.3550 236.3927 378.3173 198.8276 415.8824
## Jun 2017       311.0643 240.1020 382.0266 202.5368 419.5918
## Jul 2017       314.7736 243.8113 385.7359 206.2461 423.3011
## Aug 2017       318.4829 247.5205 389.4453 209.9553 427.0105

Htrain = H[1:48]
Htest = H[49:58]
Harma_fit <- auto.arima(Htrain)
Harma_forecast <- forecast(Harma_fit, h = 5)
Harma_fit_accuracy <- accuracy(Harma_forecast, Htest)
Harma_fit; Harma_forecast; Harma_fit_accuracy

## Series: Htrain 
## ARIMA(0,1,1)                    
## 
## Coefficients:
##           ma1
##       -0.4830
## s.e.   0.1622
## 
## sigma^2 estimated as 3577:  log likelihood=-258.6
## AIC=521.2   AICc=521.47   BIC=524.9

##    Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## 49       332.9718 256.3254 409.6183 215.7512 450.1924
## 50       332.9718 246.6883 419.2554 201.0125 464.9311
## 51       332.9718 238.0243 427.9193 187.7621 478.1815
## 52       332.9718 230.0874 435.8563 175.6237 490.3200
## 53       332.9718 222.7204 443.2233 164.3568 501.5869

##                     ME     RMSE      MAE        MPE     MAPE      MASE
## Training set  10.28526 58.54828 47.12534  -8.812042 37.80918 0.8870208
## Test set     -77.97183 86.23460 77.97183 -33.680759 33.68076        NA
##                    ACF1
## Training set 0.02400877
## Test set             NA

plot(Harma_forecast, ylim=c(0,500))
lines(H[1:58])

Hfit1 <- meanf(H,h=3)
Hfit2 <- rwf(H,h=3)
Hfit3 <- snaive(H,h=3)

plot(Hfit1, plot.conf=FALSE,
  main="Forecasts for 1/4 SOLARBAN 60 100 X 144")
lines(Hfit2$mean,col=2)
lines(Hfit3$mean,col=3)
lines(H)
legend("topleft", lty=1, col=c(4,2,3),
  legend=c("Mean method","Naive method","Seasonal naive method"))

cycle(H)

##      Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
## 2012   1   2   3   4   5   6   7   8   9  10  11  12
## 2013   1   2   3   4   5   6   7   8   9  10  11  12
## 2014   1   2   3   4   5   6   7   8   9  10  11  12
## 2015   1   2   3   4   5   6   7   8   9  10  11  12
## 2016   1   2   3   4   5   6   7   8   9  10

boxplot(H~cycle(H))

5/32 CLEAR 96X130

 I <- ts(c(157,149,168,223,230,249,293,317,269,333,209,144,184,151,144,211,253,225,270,214,177,350,213,164,112,115,100,216,226,274,270,260,276,302,216,177,118,101,153,181,233,273,318,298,298,356,298,197,145,143,217,208,243,311,349,306,251,236
), 
        start = c(2012, 1), frequency = 12)

plot(I)

seasonplot(I,ylab="Demand", xlab="Month",
           main="Seasonal plot: 5/32 CLEAR 96X130 Demand",
           year.labels=TRUE, year.labels.left =TRUE, col=1:20, pch=19)

monthplot(I,ylab="Demand",xlab="Month",xaxt="n",
          main="Seasonal deviation plot: Demand for 5/32 CLEAR 96X130")

ggseasonplot(I, year.labels=TRUE, continuous=TRUE)

ets(I)

## ETS(M,N,A) 
## 
## Call:
##  ets(y = I) 
## 
##   Smoothing parameters:
##     alpha = 0.1761 
##     gamma = 1e-04 
## 
##   Initial states:
##     l = 232.0982 
##     s=-49.4809 12.0742 100.7175 27.2913 51.4044 68.476
##            38.8245 9.8169 -15.3121 -68.2627 -95.186 -80.363
## 
##   sigma:  0.1287
## 
##      AIC     AICc      BIC 
## 651.3483 662.7769 682.2550

Ifit1 <- ets(I)
Ifit2 <- ets(I,model= "MNA")
Ideviance <- 2*c(logLik(Ifit1) - logLik(Ifit2))
Idf <- attributes(logLik(Ifit1))$Idf - attributes(logLik(Ifit2))$Idf 
#P value
1-pchisq(Ideviance,Idf)

## numeric(0)

ses(I)

##          Point Forecast     Lo 80    Hi 80       Lo 95    Hi 95
## Nov 2016       236.3287 164.92658 307.7308  127.128580 345.5288
## Dec 2016       236.3287 136.37396 336.2835   83.461116 389.1963
## Jan 2017       236.3287 114.33019 358.3272   49.748073 422.9093
## Feb 2017       236.3287  95.70038 376.9570   21.256241 451.4012
## Mar 2017       236.3287  79.26498 393.3924   -3.879532 476.5369
## Apr 2017       236.3287  64.39354 408.2639  -26.623438 499.2809
## May 2017       236.3287  50.70977 421.9476  -47.550956 520.2084
## Jun 2017       236.3287  37.96773 434.6897  -67.038226 539.6956
## Jul 2017       236.3287  25.99620 446.6612  -85.347105 558.0045
## Aug 2017       236.3287  14.67030 457.9871 -102.668580 575.3260

holt(I)

##          Point Forecast     Lo 80    Hi 80      Lo 95    Hi 95
## Nov 2016       237.9794 166.55642 309.4023 128.747409 347.2113
## Dec 2016       239.5727 139.64215 339.5033  86.742099 392.4034
## Jan 2017       241.1661 119.21528 363.1169  54.658408 427.6738
## Feb 2017       242.7595 102.19317 383.3258  27.781864 457.7371
## Mar 2017       244.3529  87.35996 401.3458   4.252947 484.4528
## Apr 2017       245.9462  74.08685 417.8056 -16.890016 508.7825
## May 2017       247.5396  61.99843 433.0808 -36.221131 531.3003
## Jun 2017       249.1330  50.84934 447.4166 -54.115688 552.3817
## Jul 2017       250.7264  40.46870 460.9840 -70.834978 572.2877
## Aug 2017       252.3197  30.73192 473.9075 -86.569579 591.2091

Itrain = I[1:48]
Itest = I[49:58]
Iarma_fit <- auto.arima(Itrain)
Iarma_forecast <- forecast(Iarma_fit, h = 5)
Iarma_fit_accuracy <- accuracy(Iarma_forecast, Itest)
Iarma_fit; Iarma_forecast; Iarma_fit_accuracy

## Series: Itrain 
## ARIMA(1,0,0) with non-zero mean 
## 
## Coefficients:
##          ar1  intercept
##       0.6313   219.1812
## s.e.  0.1098    19.8703
## 
## sigma^2 estimated as 2868:  log likelihood=-258.41
## AIC=522.83   AICc=523.37   BIC=528.44

##    Point Forecast    Lo 80    Hi 80     Lo 95    Hi 95
## 49       205.1792 136.5483 273.8100 100.21735 310.1410
## 50       210.3423 129.1811 291.5035  86.21698 334.4677
## 51       213.6016 127.9565 299.2467  82.61877 344.5844
## 52       215.6591 128.2913 303.0268  82.04163 349.2765
## 53       216.9578 128.9130 305.0026  82.30493 351.6107

##                      ME     RMSE      MAE        MPE     MAPE      MASE
## Training set   1.107556 52.42537 43.81063  -6.020614 22.02724 0.9595059
## Test set     -21.147992 42.20171 32.92422 -15.998922 20.91213        NA
##                    ACF1
## Training set 0.07182829
## Test set             NA

plot(Iarma_forecast, ylim=c(0,500))
lines(I[1:58])

Ifit1 <- meanf(I,h=3)
Ifit2 <- rwf(I,h=3)
Ifit3 <- snaive(I,h=3)

plot(Ifit1, plot.conf=FALSE,
  main="Forecasts for ")
lines(Ifit2$mean,col=2)
lines(Ifit3$mean,col=3)
lines(I)
legend("topleft", lty=1, col=c(4,2,3),
  legend=c("Mean method","Naive method","Seasonal naive method"))

cycle(I)

##      Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
## 2012   1   2   3   4   5   6   7   8   9  10  11  12
## 2013   1   2   3   4   5   6   7   8   9  10  11  12
## 2014   1   2   3   4   5   6   7   8   9  10  11  12
## 2015   1   2   3   4   5   6   7   8   9  10  11  12
## 2016   1   2   3   4   5   6   7   8   9  10

boxplot(I~cycle(I))

1/4 ENEG ADV LOW-E 96X130

 J <- ts(c(15,1,8,4,4,20,32,2,16,3,5,3,3,1,1,10,3,15,9,7,13,90,229,205,171,126,128,161,197,117,297,240,197,227,266,323,105,199,176,180,255,264,244,179,271,245,162,339,225,125,309,130,192,230,236,323,168,192
), 
        start = c(2012, 1), frequency = 12)

plot(J)

seasonplot(J,ylab="Demand", xlab="Month",
           main="Seasonal plot: 5/32 CLEAR 96X130 Demand",
           year.labels=TRUE, year.labels.left =TRUE, col=1:20, pch=19)

monthplot(J,ylab="Demand",xlab="Month",xaxt="n",
          main="Seasonal deviation plot: Demand for 1/4  ENEG ADV LOW-E 96X130")

ggseasonplot(J, year.labels=TRUE, continuous=TRUE)

ets(J)

## ETS(A,N,N) 
## 
## Call:
##  ets(y = J) 
## 
##   Smoothing parameters:
##     alpha = 0.3408 
## 
##   Initial states:
##     l = 9.2794 
## 
##   sigma:  61.4965
## 
##      AIC     AICc      BIC 
## 719.3073 719.7518 725.4887

Jfit1 <- ets(J)
Jfit2 <- ets(J,model= "ANN")
Jdeviance <- 2*c(logLik(Jfit1) - logLik(Jfit2))
Jdf <- attributes(logLik(Jfit1))$Jdf - attributes(logLik(Jfit2))$Jdf 
#P value
1-pchisq(Jdeviance,Jdf)

## numeric(0)

ses(J)

##          Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## Nov 2016       213.6017 134.7908 292.4126 93.07083 334.1325
## Dec 2016       213.6017 130.3377 296.8656 86.26047 340.9429
## Jan 2017       213.6017 126.1110 301.0923 79.79630 347.4071
## Feb 2017       213.6017 122.0793 305.1240 73.63033 353.5730
## Mar 2017       213.6017 118.2179 308.9855 67.72476 359.4786
## Apr 2017       213.6017 114.5068 312.6966 62.04914 365.1542
## May 2017       213.6017 110.9298 316.2736 56.57853 370.6248
## Jun 2017       213.6017 107.4732 319.7301 51.29220 375.9112
## Jul 2017       213.6017 104.1258 323.0776 46.17270 381.0307
## Aug 2017       213.6017 100.8777 326.3257 41.20515 385.9982

holt(J)

##          Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## Nov 2016       229.7428 152.3921 307.0936 111.4450 348.0406
## Dec 2016       233.6899 153.3201 314.0597 110.7749 356.6049
## Jan 2017       237.6369 154.3555 320.9184 110.2689 365.0049
## Feb 2017       241.5840 155.4873 327.6807 109.9105 373.2575
## Mar 2017       245.5310 156.7064 334.3556 109.6855 381.3765
## Apr 2017       249.4781 158.0050 340.9511 109.5821 389.3741
## May 2017       253.4251 159.3764 347.4739 109.5900 397.2603
## Jun 2017       257.3722 160.8146 353.9297 109.7002 405.0442
## Jul 2017       261.3192 162.3148 360.3237 109.9050 412.7335
## Aug 2017       265.2663 163.8723 366.6603 110.1975 420.3351

Jtrain = J[1:48]
Jtest = J[48:58]
Jarma_fit <- auto.arima(Jtrain)
Jarma_forecast <- forecast(Jarma_fit, h = 5)
Jarma_fit_accuracy <- accuracy(Jarma_forecast, Jtest)
Jarma_fit; Jarma_forecast; Jarma_fit_accuracy

## Series: Jtrain 
## ARIMA(0,1,1) with drift         
## 
## Coefficients:
##           ma1   drift
##       -0.6533  5.7065
## s.e.   0.1310  2.9563
## 
## sigma^2 estimated as 3282:  log likelihood=-256.21
## AIC=518.42   AICc=518.97   BIC=523.97

##    Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## 49       272.4161 198.9962 345.8359 160.1301 384.7020
## 50       278.1225 200.4141 355.8309 159.2778 396.9672
## 51       283.8290 202.0566 365.6013 158.7690 408.8890
## 52       289.5354 203.8918 375.1791 158.5548 420.5161
## 53       295.2419 205.8945 384.5893 158.5969 431.8869

##                       ME      RMSE      MAE        MPE      MAPE      MASE
## Training set  -0.6007834  55.47062 37.81405 -250.86527 268.31438 0.9256565
## Test set     -58.2289792 109.69578 92.64838  -50.36837  60.74457        NA
##                    ACF1
## Training set 0.05093936
## Test set             NA

plot(Jarma_forecast, ylim=c(0,500))
lines(J[1:58])

Jfit1 <- meanf(J,h=3)
Jfit2 <- rwf(J,h=3)
Jfit3 <- snaive(J,h=3)

plot(Jfit1, plot.conf=FALSE,
  main="Forecasts for 1/4  ENEG ADV LOW-E 96X130")
lines(Jfit2$mean,col=2)
lines(Jfit3$mean,col=3)
lines(J)
legend("topleft", lty=1, col=c(4,2,3),
  legend=c("Mean method","Naive method","Seasonal naive method"))

cycle(J)

##      Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
## 2012   1   2   3   4   5   6   7   8   9  10  11  12
## 2013   1   2   3   4   5   6   7   8   9  10  11  12
## 2014   1   2   3   4   5   6   7   8   9  10  11  12
## 2015   1   2   3   4   5   6   7   8   9  10  11  12
## 2016   1   2   3   4   5   6   7   8   9  10

boxplot(J~cycle(J))

1/4 CLEAR MIRROR 96 x 130

 K <- ts(c(51,118,171,200,160,137,155,142,180,181,196,216,235,239,205,210,171,201,221,169,235,182,209,239,249,385,274,217), 
        start = c(2014, 7), frequency = 12)

plot(K)

seasonplot(K,ylab="Demand", xlab="Month",
           main="Seasonal plot: 5/32 CLEAR 96X130 Demand",
           year.labels=TRUE, year.labels.left=TRUE, col=1:20, pch=19)

monthplot(K,ylab="Demand",xlab="Month",xaxt="n",
          main="Seasonal deviation plot: Demand for 1/4  ENEG ADV LOW-E 96X130")

ggseasonplot(K, year.labels=TRUE, continuous=TRUE)

ets(K)

## ETS(M,A,N) 
## 
## Call:
##  ets(y = K) 
## 
##   Smoothing parameters:
##     alpha = 9e-04 
##     beta  = 1e-04 
## 
##   Initial states:
##     l = 126.5158 
##     b = 4.9157 
## 
##   sigma:  0.2093
## 
##      AIC     AICc      BIC 
## 310.6141 313.3414 317.2751

Kfit1 <- ets(K)
Kfit2 <- ets(K,model= "MAN")
Kdeviance <- 2*c(logLik(Kfit1) - logLik(Kfit2))
Kdf <- attributes(logLik(Kfit1))$Kdf - attributes(logLik(Kfit2))$Kdf 
#P value
1-pchisq(Kdeviance,Kdf)

## numeric(0)

ses(K)

##          Point Forecast    Lo 80    Hi 80     Lo 95    Hi 95
## Nov 2016       249.9291 191.2063 308.6519 160.12028 339.7378
## Dec 2016       249.9291 182.1545 317.7036 146.27683 353.5813
## Jan 2017       249.9291 174.1768 325.6813 134.07592 365.7822
## Feb 2017       249.9291 166.9626 332.8955 123.04283 376.8153
## Mar 2017       249.9291 160.3274 339.5307 112.89519 386.9629
## Apr 2017       249.9291 154.1508 345.7073 103.44886 396.4092
## May 2017       249.9291 148.3491 351.5090  94.57586 405.2822
## Jun 2017       249.9291 142.8612 356.9969  86.18297 413.6751
## Jul 2017       249.9291 137.6413 362.2168  78.19977 421.6583
## Aug 2017       249.9291 132.6535 367.2046  70.57156 429.2865

holt(K)

##          Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## Nov 2016       279.3175 227.1418 331.4933 199.5216 359.1134
## Dec 2016       285.1726 232.9969 337.3484 205.3767 364.9685
## Jan 2017       291.0278 238.8520 343.2035 211.2319 370.8237
## Feb 2017       296.8829 244.7071 349.0586 217.0870 376.6788
## Mar 2017       302.7380 250.5622 354.9138 222.9421 382.5339
## Apr 2017       308.5931 256.4174 360.7689 228.7972 388.3891
## May 2017       314.4482 262.2725 366.6240 234.6523 394.2442
## Jun 2017       320.3034 268.1276 372.4792 240.5074 400.0994
## Jul 2017       326.1585 273.9827 378.3343 246.3624 405.9546
## Aug 2017       332.0136 279.8377 384.1895 252.2175 411.8097

Ktrain = K[1:21]
Ktest = K[22:29]
Karma_fit <- auto.arima(Ktrain)
Karma_forecast <- forecast(Karma_fit, h = 5)
Karma_fit_accuracy <- accuracy(Karma_forecast, Ktest)
Karma_fit; Karma_forecast; Karma_fit_accuracy

## Series: Ktrain 
## ARIMA(0,1,0)                    
## 
## sigma^2 estimated as 1214:  log likelihood=-99.39
## AIC=200.78   AICc=201.01   BIC=201.78

##    Point Forecast    Lo 80    Hi 80     Lo 95    Hi 95
## 22            235 190.3567 279.6433 166.72401 303.2760
## 23            235 171.8649 298.1351 138.44317 331.5568
## 24            235 157.6756 312.3244 116.74252 353.2575
## 25            235 145.7135 324.2865  98.44802 371.5520
## 26            235 135.1746 334.8254  82.33025 387.6698

##                     ME    RMSE      MAE       MPE     MAPE      MASE
## Training set  8.764333 33.9958 27.90719 4.5551296 16.08893 0.9524638
## Test set     17.800000 72.3837 49.40000 0.9392197 17.56365        NA
##                     ACF1
## Training set -0.09425765
## Test set              NA

plot(Karma_forecast, ylim=c(0,500))
lines(K[1:29])

Kfit1 <- meanf(K,h=3)
Kfit2 <- rwf(K,h=3)
Kfit3 <- snaive(K,h=3)

plot(Kfit1, plot.conf=FALSE,
  main="Forecasts for 1/4  ENEG ADV LOW-E 96X130")
lines(Kfit2$mean,col=2)
lines(Kfit3$mean,col=3)
lines(K)
legend("topleft", lty=1, col=c(4,2,3),
  legend=c("Mean method","Naive method","Seasonal naive method"))

cycle(K)

##      Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
## 2014                           7   8   9  10  11  12
## 2015   1   2   3   4   5   6   7   8   9  10  11  12
## 2016   1   2   3   4   5   6   7   8   9  10

boxplot(K~cycle(K))

L <- ts(c(64,186,176,185,186,191,141,72,92,62,64,139,208,162,188,198,194,217,124,89,71,84,53,140,203,214,221,226,262,239,115,100,64,45,97,195,192,257,235,250,275,256,208,134,102,82,173,146,232,263,215,262,261,248
), 
        start = c(2012, 5), frequency = 12)

plot(L)

seasonplot(L,ylab="Demand", xlab="Month",
           main="Seasonal plot: 3/32  ENEG ADV LOW-E 72X84",
           year.labels=TRUE, year.labels.left =TRUE, col=1:20, pch=19)

monthplot(L,ylab="Demand",xlab="Month",xaxt="n",
          main="Seasonal deviation plot: Demand for 3/32  ENEG ADV LOW-E 72X84")

ggseasonplot(L, year.labels=TRUE, continuous=TRUE)

ets(L)

## ETS(A,A,A) 
## 
## Call:
##  ets(y = L) 
## 
##   Smoothing parameters:
##     alpha = 1e-04 
##     beta  = 1e-04 
##     gamma = 1e-04 
## 
##   Initial states:
##     l = 114.0429 
##     b = 1.6887 
##     s=-13.6711 -70.6222 -97.4467 -82.0822 -63.081 -13.6397
##            69.1101 77.0645 59.7065 50.1901 55.6766 28.7951
## 
##   sigma:  23.6311
## 
##      AIC     AICc      BIC 
## 590.9618 607.9618 624.7746

Lfit1 <- ets(L)
Lfit2 <- ets(L,model= "MNN")
Ldeviance <- 2*c(logLik(Lfit1) - logLik(Lfit2))
Ldf <- attributes(logLik(Lfit1))$Ldf - attributes(logLik(Lfit2))$Ldf 
#P value
1-pchisq(Ldeviance,Ldf)

## numeric(0)

ses(L)

##          Point Forecast     Lo 80    Hi 80      Lo 95    Hi 95
## Nov 2016       248.0013 185.33558 310.6670 152.162345 343.8403
## Dec 2016       248.0013 159.38301 336.6196 112.471328 383.5313
## Jan 2017       248.0013 139.46832 356.5343  82.014429 413.9882
## Feb 2017       248.0013 122.67925 373.3233  56.337768 439.6648
## Mar 2017       248.0013 107.88769 388.1149  33.716029 462.2866
## Apr 2017       248.0013  94.51504 401.4876  13.264329 482.7383
## May 2017       248.0013  82.21759 413.7850  -5.543003 501.5456
## Jun 2017       248.0013  70.77138 425.2312 -23.048477 519.0511
## Jul 2017       248.0013  60.02084 435.9818 -39.490003 535.4926
## Aug 2017       248.0013  49.85272 446.1499 -55.040805 551.0434

holt(L)

##          Point Forecast    Lo 80    Hi 80      Lo 95    Hi 95
## Nov 2016       256.0476 191.1233 320.9718 156.754453 355.3407
## Dec 2016       264.0930 171.8744 356.3116 123.056769 405.1292
## Jan 2017       272.1384 158.8643 385.4126  98.900549 445.3763
## Feb 2017       280.1839 149.0519 411.3159  79.634754 480.7330
## Mar 2017       288.2293 141.2662 435.1924  63.468652 512.9900
## Apr 2017       296.2747 134.9092 457.6403  49.487464 543.0620
## May 2017       304.3202 129.6267 479.0137  37.149521 571.4908
## Jun 2017       312.3656 125.1885 499.5427  26.102913 598.6283
## Jul 2017       320.4111 121.4353 519.3868  16.103857 624.7183
## Aug 2017       328.4565 118.2512 538.6617   6.975311 649.9377

Ltrain = L[1:43]
Ltest = L[44:53]
Larma_fit <- auto.arima(Ltrain)
Larma_forecast <- forecast(Larma_fit, h = 5)
Larma_fit_accuracy <- accuracy(Larma_forecast, Ltest)
Larma_fit; Larma_forecast; Larma_fit_accuracy

## Series: Ltrain 
## ARIMA(2,0,1) with non-zero mean 
## 
## Coefficients:
##          ar1      ar2      ma1  intercept
##       1.5867  -0.8101  -0.7242   154.4992
## s.e.  0.1094   0.0993   0.1293     8.1128
## 
## sigma^2 estimated as 1750:  log likelihood=-220.35
## AIC=450.7   AICc=452.32   BIC=459.51

##    Point Forecast     Lo 80    Hi 80     Lo 95    Hi 95
## 44      146.28295 92.665670 199.9002  64.28240 228.2835
## 45       98.11922 27.310497 168.9279 -10.17338 206.4118
## 46       71.69473 -5.187417 148.5769 -45.88637 189.2758
## 47       68.78490 -8.751850 146.3216 -49.79732 187.3671
## 48       85.57513  7.594319 163.5560 -33.68623 204.8365

##                      ME     RMSE      MAE       MPE     MAPE      MASE
## Training set -0.1772618 39.84435 32.83674 -9.938478 27.14797 0.8990502
## Test set     33.3086140 54.37662 38.22179 21.766509 25.43306        NA
##                    ACF1
## Training set -0.1486797
## Test set             NA

plot(Larma_forecast, ylim=c(0,500))
lines(L[1:58])

Lfit1 <- meanf(L,h=3)
Lfit2 <- rwf(L,h=3)
Lfit3 <- snaive(L,h=3)

plot(Lfit1, plot.conf=FALSE,
  main="Forecasts for 3/32  ENEG ADV LOW-E 72X84")
lines(Lfit2$mean,col=2)
lines(Lfit3$mean,col=3)
lines(L)
legend("topleft", lty=1, col=c(4,2,3),
  legend=c("Mean method","Naive method","Seasonal naive method"))

cycle(L)

##      Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
## 2012                   5   6   7   8   9  10  11  12
## 2013   1   2   3   4   5   6   7   8   9  10  11  12
## 2014   1   2   3   4   5   6   7   8   9  10  11  12
## 2015   1   2   3   4   5   6   7   8   9  10  11  12
## 2016   1   2   3   4   5   6   7   8   9  10

boxplot(L~cycle(L))

 M <- ts(c(73,79,102,112,122,114,96,125,83,98,80,96,96,101,97,147,145,154,137,122,101,106,92,108,94,133,95,124,145,162,137,123,146,135,96,166,119,125,153,184,175,204,213,150,154,134,148,58,19,29,29,27,42,155,143,184,210,180
), 
        start = c(2012, 1), frequency = 12)

plot(M)

seasonplot(M,ylab="Demand", xlab="Month",
           main="Seasonal plot: 3/8 SHOWERGUARD 96X130 Demand",
           year.labels=TRUE, year.labels.left =TRUE, col=1:20, pch=19)

monthplot(M,ylab="Demand",xlab="Month",xaxt="n",
          main="Seasonal deviation plot: Demand for 3/8 SHOWERGUARD 96X130")

ggseasonplot(M, year.labels=TRUE, continuous=TRUE)

ets(M)

## ETS(A,N,N) 
## 
## Call:
##  ets(y = M) 
## 
##   Smoothing parameters:
##     alpha = 0.8564 
## 
##   Initial states:
##     l = 74.3316 
## 
##   sigma:  30.9283
## 
##      AIC     AICc      BIC 
## 639.5796 640.0241 645.7610

Mfit1 <- ets(M)
Mfit2 <- ets(M,model= "MNN")
Mdeviance <- 2*c(logLik(Mfit1) - logLik(Mfit2))
Mdf <- attributes(logLik(Mfit1))$Mdf - attributes(logLik(Mfit2))$Mdf 
#P value
1-pchisq(Mdeviance,Mdf)

## numeric(0)

ses(M)

##          Point Forecast     Lo 80    Hi 80     Lo 95    Hi 95
## Nov 2016       183.6414 144.00522 223.2777 123.02308 244.2598
## Dec 2016       183.6414 131.44745 235.8354 103.81762 263.4653
## Jan 2017       183.6414 121.37270 245.9102  88.40962 278.8733
## Feb 2017       183.6414 112.71486 254.5680  75.16860 292.1143
## Mar 2017       183.6414 105.00454 262.2783  63.37668 303.9062
## Apr 2017       183.6414  97.98547 269.2974  52.64194 314.6409
## May 2017       183.6414  91.49955 275.7833  42.72257 324.5603
## Jun 2017       183.6414  85.44108 281.8418  33.45694 333.8259
## Jul 2017       183.6414  79.73526 287.5476  24.73065 342.5522
## Aug 2017       183.6414  74.32686 292.9560  16.45921 350.8237

holt(M)

##          Point Forecast     Lo 80    Hi 80     Lo 95    Hi 95
## Nov 2016       186.0171 146.36683 225.6674 125.37725 246.6570
## Dec 2016       187.9084 135.88540 239.9315 108.34607 267.4708
## Jan 2017       189.7998 127.82465 251.7749  95.01701 284.5825
## Feb 2017       191.6911 121.15245 262.2297  83.81155 299.5706
## Mar 2017       193.5824 115.41111 271.7537  74.02972 313.1351
## Apr 2017       195.4737 110.34988 280.5976  65.28802 325.6595
## May 2017       197.3651 105.81367 288.9165  57.34928 337.3809
## Jun 2017       199.2564 101.69865 296.8142  50.05470 348.4581
## Jul 2017       201.1477  97.93128 304.3642  43.29179 359.0037
## Aug 2017       203.0391  94.45718 311.6210  36.97740 369.1007

Mtrain = M[1:48]
Mtest = M[49:58]
Marma_fit <- auto.arima(Mtrain)
Marma_forecast <- forecast(Marma_fit, h = 5)
Marma_fit_accuracy <- accuracy(Marma_forecast, Mtest)
Marma_fit; Marma_forecast; Marma_fit_accuracy

## Series: Mtrain 
## ARIMA(1,1,0)                    
## 
## Coefficients:
##           ar1
##       -0.4206
## s.e.   0.1486
## 
## sigma^2 estimated as 721.1:  log likelihood=-220.93
## AIC=445.86   AICc=446.13   BIC=449.56

##    Point Forecast    Lo 80    Hi 80     Lo 95    Hi 95
## 49       95.85701 61.44390 130.2701 43.226694 148.4873
## 50       79.93309 40.16150 119.7047 19.107701 140.7585
## 51       86.63122 39.10057 134.1619 13.939362 159.3231
## 52       83.81376 30.80631 136.8212  2.745853 164.8817
## 53       84.99888 26.58538 143.4124 -4.336866 174.3346

##                       ME     RMSE      MAE         MPE     MAPE      MASE
## Training set   0.3346652 26.28731 20.38758   -3.145108  17.7117 0.9187115
## Test set     -57.0467931 58.13647 57.04679 -218.333982 218.3340        NA
##                      ACF1
## Training set -0.004227719
## Test set               NA

plot(Marma_forecast, ylim=c(0,500))
lines(M[1:58])

Mfit1 <- meanf(M,h=3)
Mfit2 <- rwf(M,h=3)
Mfit3 <- snaive(M,h=3)

plot(Mfit1, plot.conf=FALSE,
  main="Forecasts for 3/8 SHOWERGUARD 96X130")
lines(Mfit2$mean,col=2)
lines(Mfit3$mean,col=3)
lines(M)
legend("topleft", lty=1, col=c(4,2,3),
  legend=c("Mean method","Naive method","Seasonal naive method"))

cycle(M)

##      Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
## 2012   1   2   3   4   5   6   7   8   9  10  11  12
## 2013   1   2   3   4   5   6   7   8   9  10  11  12
## 2014   1   2   3   4   5   6   7   8   9  10  11  12
## 2015   1   2   3   4   5   6   7   8   9  10  11  12
## 2016   1   2   3   4   5   6   7   8   9  10

boxplot(M~cycle(M))

 N <- ts(c(2,26,7,6,3,2,6,9,37,19,11,5,5,15,20,27,26,29,34,54,39,46,38,23,6,24,11,45,35,45,40,79,52,62,30,25,48,52,1,78,70,110,91,125,176,155,130,167,95,19,126,137,150,158,150,210,233,148
), 
        start = c(2012, 1), frequency = 12)

plot(N)

seasonplot(N,ylab="Demand", xlab="Month",
           main="Seasonal plot: 1/8 SOLARBAN60 TC 72X84 Demand",
           year.labels=TRUE, year.labels.left =TRUE, col=1:20, pch=19)

monthplot(N,ylab="Demand",xlab="Month",xaxt="n",
          main="Seasonal deviation plot: Demand for 1/8 SOLARBAN60 TC 72X84")

ggseasonplot(N, year.labels=TRUE, continuous=TRUE)

ets(N)

## ETS(M,A,N) 
## 
## Call:
##  ets(y = N) 
## 
##   Smoothing parameters:
##     alpha = 0.1787 
##     beta  = 1e-04 
## 
##   Initial states:
##     l = 13.365 
##     b = 1.398 
## 
##   sigma:  0.5653
## 
##      AIC     AICc      BIC 
## 612.7938 613.9477 623.0960

Nfit1 <- ets(N)
Nfit2 <- ets(N,model= "MAN")
Ndeviance <- 2*c(logLik(Nfit1) - logLik(Nfit2))
Ndf <- attributes(logLik(Nfit1))$Ndf - attributes(logLik(Nfit2))$Ndf 
#P value
1-pchisq(Ndeviance,Ndf)

## numeric(0)

ses(N)

##          Point Forecast     Lo 80    Hi 80     Lo 95    Hi 95
## Nov 2016       175.1096 136.14063 214.0787 115.51168 234.7076
## Dec 2016       175.1096 130.01670 220.2026 106.14594 244.0733
## Jan 2017       175.1096 124.63032 225.5890  97.90817 252.3111
## Feb 2017       175.1096 119.76570 230.4536  90.46839 259.7509
## Mar 2017       175.1096 115.29542 234.9239  83.63169 266.5876
## Apr 2017       175.1096 111.13676 239.0825  77.27156 272.9477
## May 2017       175.1096 107.23240 242.9869  71.30037 278.9189
## Jun 2017       175.1096 103.54073 246.6785  65.65444 284.5648
## Jul 2017       175.1096 100.03036 250.1889  60.28580 289.9335
## Aug 2017       175.1096  96.67695 253.5423  55.15719 295.0621

holt(N)

##          Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## Nov 2016       182.9910 144.6390 221.3430 124.3367 241.6453
## Dec 2016       186.0416 143.0821 229.0011 120.3408 251.7425
## Jan 2017       189.0922 141.9721 236.2123 117.0282 261.1562
## Feb 2017       192.1428 141.1992 243.0864 114.2313 270.0543
## Mar 2017       195.1934 140.6926 249.6942 111.8417 278.5452
## Apr 2017       198.2440 140.4031 256.0849 109.7840 286.7041
## May 2017       201.2946 140.2950 262.2943 108.0037 294.5855
## Jun 2017       204.3452 140.3414 268.3490 106.4598 302.2306
## Jul 2017       207.3958 140.5215 274.2701 105.1204 309.6712
## Aug 2017       210.4464 140.8188 280.0740 103.9602 316.9327

Ntrain = N[1:48]
Ntest = N[49:58]
Narma_fit <- auto.arima(Ntrain)
Narma_forecast <- forecast(Narma_fit, h = 5)
Narma_fit_accuracy <- accuracy(Narma_forecast, Ntest)
Narma_fit; Narma_forecast; Narma_fit_accuracy

## Series: Ntrain 
## ARIMA(1,1,0)                    
## 
## Coefficients:
##           ar1
##       -0.4000
## s.e.   0.1372
## 
## sigma^2 estimated as 481.4:  log likelihood=-211.42
## AIC=426.85   AICc=427.12   BIC=430.55

##    Point Forecast    Lo 80    Hi 80     Lo 95    Hi 95
## 49       152.2006 124.0819 180.3192 109.19683 195.2043
## 50       158.1201 125.3282 190.9120 107.96922 208.2710
## 51       155.7524 116.6117 194.8931  95.89180 215.6130
## 52       156.6994 112.9385 200.4603  89.77293 223.6260
## 53       156.3206 108.0845 204.5568  82.54984 230.0914

##                      ME     RMSE      MAE       MPE     MAPE      MASE
## Training set   4.462428 21.47912 15.37082 -117.7867 157.7703 0.8799374
## Test set     -50.418631 69.19470 50.41863 -166.9256 166.9256        NA
##                     ACF1
## Training set -0.05600483
## Test set              NA

plot(Narma_forecast, ylim=c(0,500))
lines(N[1:58])

Nfit1 <- meanf(N,h=3)
Nfit2 <- rwf(N,h=3)
Nfit3 <- snaive(N,h=3)

plot(Nfit1, plot.conf=FALSE,
  main="Forecasts for 1/8 SOLARBAN60 TC 72X84")
lines(Nfit2$mean,col=2)
lines(Nfit3$mean,col=3)
lines(N)
legend("topleft", lty=1, col=c(4,2,3),
  legend=c("Mean method","Naive method","Seasonal naive method"))

cycle(N)

##      Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
## 2012   1   2   3   4   5   6   7   8   9  10  11  12
## 2013   1   2   3   4   5   6   7   8   9  10  11  12
## 2014   1   2   3   4   5   6   7   8   9  10  11  12
## 2015   1   2   3   4   5   6   7   8   9  10  11  12
## 2016   1   2   3   4   5   6   7   8   9  10

boxplot(N~cycle(N))

 O <- ts(c(3,8,5,3,22,2,4,3,2,25,2,19,1,7,2,18,12,7,62,43,64,92,57,45,27,70,44,72,78,104,110,73,111,66,51,50,51,28,55,76,96,135,75,92,118,247,156,94,96,71,141,86,124,96,107,104,113,170
), 
        start = c(2012, 1), frequency = 12)

plot(O)

seasonplot(O,ylab="Demand", xlab="Month",
           main="Seasonal plot: 1/8 SOLARBAN60 72X84 Demand",
           year.labels=TRUE, year.labels.left =TRUE, col=1:20, pch=19)

monthplot(O,ylab="Demand",xlab="Month",xaxt="n",
          main="Seasonal deviation plot: Demand for 1/8 SOLARBAN60 72X84")

ggseasonplot(O, year.labels=TRUE, continuous=TRUE)

ets(O)

## ETS(M,A,N) 
## 
## Call:
##  ets(y = O) 
## 
##   Smoothing parameters:
##     alpha = 0.04 
##     beta  = 1e-04 
## 
##   Initial states:
##     l = 1.1744 
##     b = 1.3739 
## 
##   sigma:  0.7416
## 
##      AIC     AICc      BIC 
## 629.7916 630.9454 640.0938

Ofit1 <- ets(O)
Ofit2 <- ets(O,model= "MNN")
Odeviance <- 2*c(logLik(Ofit1) - logLik(Ofit2))
Odf <- attributes(logLik(Ofit1))$Odf - attributes(logLik(Ofit2))$Odf 
#P value
1-pchisq(Odeviance,Odf)

## numeric(0)

ses(O)

##          Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## Nov 2016       138.0484 96.99953 179.0973 75.26955 200.8273
## Dec 2016       138.0484 92.54854 183.5483 68.46235 207.6345
## Jan 2017       138.0484 88.49575 187.6011 62.26414 213.8327
## Feb 2017       138.0484 84.75025 191.3466 56.53590 219.5610
## Mar 2017       138.0484 81.25122 194.8457 51.18459 224.9123
## Apr 2017       138.0484 77.95558 198.1413 46.14434 229.9525
## May 2017       138.0484 74.83151 201.2654 41.36649 234.7304
## Jun 2017       138.0484 71.85473 204.2421 36.81389 239.2830
## Jul 2017       138.0484 69.00617 207.0907 32.45740 243.6395
## Aug 2017       138.0484 66.27056 209.8263 28.27366 247.8232

holt(O)

##          Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## Nov 2016       138.0276 97.53300 178.5221 76.09648 199.9586
## Dec 2016       140.3502 96.97516 183.7253 74.01378 206.6867
## Jan 2017       142.6729 96.59568 188.7502 72.20386 213.1420
## Feb 2017       144.9956 96.36482 193.6264 70.62124 219.3700
## Mar 2017       147.3183 96.26029 198.3763 69.23181 225.4048
## Apr 2017       149.6410 96.26483 203.0171 68.00921 231.2728
## May 2017       151.9637 96.36482 207.5625 66.93257 236.9948
## Jun 2017       154.2864 96.54922 212.0235 65.98502 242.5877
## Jul 2017       156.6090 96.80897 216.4091 65.15273 248.0653
## Aug 2017       158.9317 97.13654 220.7269 64.42414 253.4393

Otrain = O[1:48]
Otest = O[49:58]
Oarma_fit <- auto.arima(Otrain)
Oarma_forecast <- forecast(Oarma_fit, h = 5)
Oarma_fit_accuracy <- accuracy(Oarma_forecast, Otest)
Oarma_fit; Oarma_forecast; Oarma_fit_accuracy

## Series: Otrain 
## ARIMA(0,1,1)                    
## 
## Coefficients:
##           ma1
##       -0.4930
## s.e.   0.1607
## 
## sigma^2 estimated as 1058:  log likelihood=-229.99
## AIC=463.97   AICc=464.25   BIC=467.68

##    Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## 49       129.6632 87.97263 171.3537 65.90298 193.4234
## 50       129.6632 82.92058 176.4058 58.17655 201.1498
## 51       129.6632 78.36368 180.9627 51.20737 208.1190
## 52       129.6632 74.17979 185.1466 44.80866 214.5177
## 53       129.6632 70.29000 189.0364 38.85974 220.4666

##                      ME     RMSE      MAE       MPE     MAPE     MASE
## Training set   5.175388 31.84636 21.71630 -59.17987 98.18645 0.889857
## Test set     -26.063175 36.44625 30.59791 -32.99758 36.21370       NA
##                    ACF1
## Training set 0.05031384
## Test set             NA

plot(Oarma_forecast, ylim=c(0,300))
lines(O[1:58])

Ofit1 <- meanf(O,h=3)
Ofit2 <- rwf(O,h=3)
Ofit3 <- snaive(O,h=3)

plot(Ofit1, plot.conf=FALSE,
  main="Forecasts for 1/8 SOLARBAN60 72X84")
lines(Ofit2$mean,col=2)
lines(Ofit3$mean,col=3)
lines(O)
legend("topleft", lty=1, col=c(4,2,3),
  legend=c("Mean method","Naive method","Seasonal naive method"))

cycle(O)

##      Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
## 2012   1   2   3   4   5   6   7   8   9  10  11  12
## 2013   1   2   3   4   5   6   7   8   9  10  11  12
## 2014   1   2   3   4   5   6   7   8   9  10  11  12
## 2015   1   2   3   4   5   6   7   8   9  10  11  12
## 2016   1   2   3   4   5   6   7   8   9  10

boxplot(O~cycle(O))

1/4 SOLARBAN 70XL TC 100x144

 P <- ts(c(3,9,50,28,45,30,26,28,115,17,5,36,13,1,76,32,92,86,50,65,37,44,79,19,99,78,96,30,50,109,41,46,97,255,175,140,134,193,85,178,50,144,68,111), 
        start = c(2013, 3), frequency = 12)

plot(P)

seasonplot(P,ylab="Demand", xlab="Month",
           main="Seasonal plot: 1/4 SOLARBAN 70XL TC 100x144 Demand",
           year.labels=TRUE, year.labels.left =TRUE, col=1:20, pch=19)

monthplot(P,ylab="Demand",xlab="Month",xaxt="n",
          main="Seasonal deviation plot: Demand for 1/4 SOLARBAN 70XL TC 100x144")

ggseasonplot(P, year.labels=TRUE, continuous=TRUE)

ets(P)

## ETS(M,Ad,N) 
## 
## Call:
##  ets(y = P) 
## 
##   Smoothing parameters:
##     alpha = 6e-04 
##     beta  = 1e-04 
##     phi   = 0.9673 
## 
##   Initial states:
##     l = 9.862 
##     b = 4.0461 
## 
##   sigma:  0.658
## 
##      AIC     AICc      BIC 
## 503.5232 505.7934 514.2283

Pfit1 <- ets(P)
Pfit2 <- ets(P,model= "ANN")
Pdeviance <- 2*c(logLik(Pfit1) - logLik(Pfit2))
Pdf <- attributes(logLik(Pfit1))$Pdf - attributes(logLik(Pfit2))$Pdf 
#P value
1-pchisq(Pdeviance,Pdf)

## numeric(0)

ses(P)

##          Point Forecast    Lo 80    Hi 80      Lo 95    Hi 95
## Nov 2016       111.1485 51.44544 170.8516 19.8405353 202.4565
## Dec 2016       111.1485 49.49309 172.8040 16.8546772 205.4424
## Jan 2017       111.1485 47.60070 174.6964 13.9605088 208.3366
## Feb 2017       111.1485 45.76305 176.5340 11.1500688 211.1470
## Mar 2017       111.1485 43.97566 178.3214  8.4164852 213.8806
## Apr 2017       111.1485 42.23461 180.0625  5.7537777 216.5433
## May 2017       111.1485 40.53647 181.7606  3.1567035 219.1404
## Jun 2017       111.1485 38.87823 183.4188  0.6206359 221.6764
## Jul 2017       111.1485 37.25718 185.0399 -1.8585325 224.1556
## Aug 2017       111.1485 35.67095 186.6261 -4.2844677 226.5815

holt(P)

##          Point Forecast     Lo 80    Hi 80    Lo 95    Hi 95
## Nov 2016       133.1835  77.80859 188.5584 48.49489 217.8720
## Dec 2016       135.8719  80.49700 191.2468 51.18331 220.5605
## Jan 2017       138.5603  83.18541 193.9352 53.87172 223.2489
## Feb 2017       141.2487  85.87383 196.6236 56.56013 225.9373
## Mar 2017       143.9371  88.56223 199.3120 59.24853 228.6258
## Apr 2017       146.6256  91.25064 202.0005 61.93692 231.3142
## May 2017       149.3140  93.93904 204.6889 64.62532 234.0026
## Jun 2017       152.0024  96.62743 207.3774 67.31370 236.6911
## Jul 2017       154.6908  99.31582 210.0658 70.00207 239.3795
## Aug 2017       157.3792 102.00421 212.7543 72.69044 242.0680

Ptrain = P[1:35]
Ptest = P[36:44]
Parma_fit <- auto.arima(Ptrain)
Parma_forecast <- forecast(Parma_fit, h = 5)
Parma_fit_accuracy <- accuracy(Parma_forecast, Ptest)
Parma_fit; Parma_forecast; Parma_fit_accuracy

## Series: Ptrain 
## ARIMA(0,1,1)                    
## 
## Coefficients:
##           ma1
##       -0.6023
## s.e.   0.1932
## 
## sigma^2 estimated as 2193:  log likelihood=-178.74
## AIC=361.48   AICc=361.87   BIC=364.53

##    Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## 36       156.9502 96.93941 216.9609 65.17163 248.7287
## 37       156.9502 92.36687 221.5335 58.17855 255.7218
## 38       156.9502 88.09734 225.8030 51.64885 262.2515
## 39       156.9502 84.07752 229.8228 45.50107 268.3993
## 40       156.9502 80.26814 233.6322 39.67513 274.2252

##                     ME     RMSE      MAE        MPE      MAPE      MASE
## Training set  10.72867 45.46906 31.24535 -100.25116 144.41495 0.7799867
## Test set     -10.95016 39.32818 33.79003  -16.67541  28.87718        NA
##                      ACF1
## Training set -0.003023493
## Test set               NA

plot(Parma_forecast, ylim=c(0,300))
lines(P[1:44])

Pfit1 <- meanf(P,h=3)
Pfit2 <- rwf(P,h=3)
Pfit3 <- snaive(P,h=3)

plot(Pfit1, plot.conf=FALSE,
  main="Forecasts for 1/4 SOLARBAN 70XL TC 100x144")
lines(Pfit2$mean,col=2)
lines(Pfit3$mean,col=3)
lines(P)
legend("topleft", lty=1, col=c(4,2,3),
  legend=c("Mean method","Naive method","Seasonal naive method"))

cycle(P)

##      Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
## 2013           3   4   5   6   7   8   9  10  11  12
## 2014   1   2   3   4   5   6   7   8   9  10  11  12
## 2015   1   2   3   4   5   6   7   8   9  10  11  12
## 2016   1   2   3   4   5   6   7   8   9  10

boxplot(P~cycle(P))

1/4 STARPHIRE 96X130

 Q <- ts(c(24,3,12,25,1,10,22,14,9,10,9,14,9,5,7,16,42,28,20,20,10,27,25,45,6,10,14,10,41,103,22,82,101,154,105,104,85,41,75,72,59,135,250,84,130,144,139,94,157,52,70,52,56,73,89,123,139,143
), 
        start = c(2012, 1), frequency = 12)

plot(Q)

seasonplot(Q,ylab="Demand", xlab="Month",
           main="Seasonal plot: 1/4 STARPHIRE 96X130 Demand",
           year.labels=TRUE, year.labels.left =TRUE, col=1:20, pch=19)

monthplot(Q,ylab="Demand",xlab="Month",xaxt="n",
          main="Seasonal deviation plot: Demand for 1/4 STARPHIRE 96X130")

ggseasonplot(Q, year.labels=TRUE, continuous=TRUE)

ets(Q)

## ETS(M,A,N) 
## 
## Call:
##  ets(y = Q) 
## 
##   Smoothing parameters:
##     alpha = 0.8647 
##     beta  = 1e-04 
## 
##   Initial states:
##     l = 12.1355 
##     b = 4.9559 
## 
##   sigma:  0.5666
## 
##      AIC     AICc      BIC 
## 616.7643 617.9181 627.0665

Qfit1 <- ets(Q)
Qfit2 <- ets(Q,model= "MNN")
Qdeviance <- 2*c(logLik(Qfit1) - logLik(Qfit2))
Qdf <- attributes(logLik(Qfit1))$Qdf - attributes(logLik(Qfit2))$Qdf 
#P value
1-pchisq(Qdeviance,Qdf)

## numeric(0)

ses(Q)

##          Point Forecast    Lo 80    Hi 80     Lo 95    Hi 95
## Nov 2016       128.7272 81.99811 175.4564 57.261202 200.1933
## Dec 2016       128.7272 77.38057 180.0739 50.199288 207.2552
## Jan 2017       128.7272 73.14533 184.3091 43.722041 213.7324
## Feb 2017       128.7272 69.21071 188.2438 37.704560 219.7499
## Mar 2017       128.7272 65.52054 191.9339 32.060943 225.3935
## Apr 2017       128.7272 62.03425 195.4202 26.729115 230.7253
## May 2017       128.7272 58.72136 198.7331 21.662484 235.7920
## Jun 2017       128.7272 55.55831 201.8962 16.825022 240.6294
## Jul 2017       128.7272 52.52645 204.9280 12.188186 245.2663
## Aug 2017       128.7272 49.61068 207.8438  7.728911 249.7256

holt(Q)

##          Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## Nov 2016       131.2286 84.85781 177.5995 60.31057 202.1467
## Dec 2016       133.2676 82.96438 183.5707 56.33549 210.1996
## Jan 2017       135.3065 81.35505 189.2579 52.79489 217.8180
## Feb 2017       137.3454 79.97558 194.7152 49.60584 225.0849
## Mar 2017       139.3843 78.78706 199.9815 46.70883 232.0597
## Apr 2017       141.4232 77.76044 205.0859 44.05942 238.7870
## May 2017       143.4621 76.87337 210.0508 41.62342 245.3008
## Jun 2017       145.5010 76.10818 214.8939 39.37383 251.6282
## Jul 2017       147.5399 75.45065 219.6292 37.28889 257.7910
## Aug 2017       149.5788 74.88911 224.2686 35.35076 263.8069

Qtrain = Q[1:48]
Qtest = Q[49:58]
Qarma_fit <- auto.arima(Qtrain)
Qarma_forecast <- forecast(Qarma_fit, h = 5)
Qarma_fit_accuracy <- accuracy(Qarma_forecast, Qtest)
Qarma_fit; Qarma_forecast; Qarma_fit_accuracy

## Series: Qtrain 
## ARIMA(0,1,1)                    
## 
## Coefficients:
##           ma1
##       -0.6234
## s.e.   0.1348
## 
## sigma^2 estimated as 1326:  log likelihood=-235.4
## AIC=474.81   AICc=475.08   BIC=478.51

##    Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## 49       120.1023 73.42686 166.7778 48.71836 191.4863
## 50       120.1023 70.22651 169.9781 43.82385 196.3808
## 51       120.1023 67.21948 172.9852 39.22500 200.9796
## 52       120.1023 64.37448 175.8302 34.87394 205.3307
## 53       120.1023 61.66782 178.5368 30.73447 209.4702

##                      ME     RMSE      MAE       MPE     MAPE      MASE
## Training set   5.576082 35.65421 20.56028 -61.80016 94.67156 0.7973045
## Test set     -42.702317 58.74800 57.46139 -84.89469 94.29537        NA
##                    ACF1
## Training set 0.02468911
## Test set             NA

plot(Qarma_forecast, ylim=c(0,300))
lines(Q[1:58])

Qfit1 <- meanf(Q,h=3)
Qfit2 <- rwf(Q,h=3)
Qfit3 <- snaive(Q,h=3)

plot(Qfit1, plot.conf=FALSE,
  main="Forecasts for 1/4 STARPHIRE 96X130")
lines(Qfit2$mean,col=2)
lines(Qfit3$mean,col=3)
lines(Q)
legend("topleft", lty=1, col=c(4,2,3),
  legend=c("Mean method","Naive method","Seasonal naive method"))

cycle(Q)

##      Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
## 2012   1   2   3   4   5   6   7   8   9  10  11  12
## 2013   1   2   3   4   5   6   7   8   9  10  11  12
## 2014   1   2   3   4   5   6   7   8   9  10  11  12
## 2015   1   2   3   4   5   6   7   8   9  10  11  12
## 2016   1   2   3   4   5   6   7   8   9  10

boxplot(Q~cycle(Q))

Twelve <- read.csv("Twelve.csv", header = TRUE, stringsAsFactors = FALSE)
Thirteen <- read.csv("Thirteen.csv", header = TRUE, stringsAsFactors = FALSE )
Fourteen <- read.csv("Fourteen.csv", header = TRUE, stringsAsFactors = FALSE)
TwelveThirteen <- inner_join(Twelve, Thirteen, by= "Item_Code")
kable(head(TwelveThirteen))

Item_Code	X1.1.2012	X2.1.2012	X3.1.2012	X4.1.2012	X5.1.2012	X6.1.2012	X7.1.2012	X8.1.2012	X9.1.2012	X10.1.2012	X11.1.2012	X12.1.2012	X1.1.2013	X2.1.2013	X3.1.2013	X4.1.2013	X5.1.2013	X6.1.2013	X7.1.2013	X8.1.2013	X9.1.2013	X10.1.2013	X11.1.2013	X12.1.2013
SS02CLR3636	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
SS02CLR3660	18	3	5	1	1	2	1	6	0	2	5	6	12	8	3	7	3	2	0	108	1	3	0	3
SS02CLR4884	2	0	4	6	2	7	0	1	0	92	0	1	3	5	1	0	0	1	1	35	157	1	0	1
SS02CLR7284	468	538	687	932	1,162	1,067	1,131	1,201	1,183	1,082	793	505	498	382	401	752	955	878	999	1,053	1,007	1,283	753	396
SS02LOE7284	0	0	0	0	64	186	176	185	186	191	141	72	92	62	64	139	208	162	188	198	194	217	124	89
SS03ACC7284	30	33	44	48	117	93	133	126	140	101	48	37	29	24	23	40	79	39	57	28	18	15	9	6

TTF <- inner_join(TwelveThirteen, Fourteen, by = "Item_Code")
kable(head(TTF))

Item_Code	X1.1.2012	X2.1.2012	X3.1.2012	X4.1.2012	X5.1.2012	X6.1.2012	X7.1.2012	X8.1.2012	X9.1.2012	X10.1.2012	X11.1.2012	X12.1.2012	X1.1.2013	X2.1.2013	X3.1.2013	X4.1.2013	X5.1.2013	X6.1.2013	X7.1.2013	X8.1.2013	X9.1.2013	X10.1.2013	X11.1.2013	X12.1.2013	X1.1.2014	X2.1.2014	X3.1.2014	X4.1.2014	X5.1.2014	X6.1.2014	X7.1.2014	X8.1.2014	X9.1.2014	X10.1.2014	X11.1.2014	X12.1.2014
SS02CLR3660	18	3	5	1	1	2	1	6	0	2	5	6	12	8	3	7	3	2	0	108	1	3	0	3	2	3	1	0	1	0	0	0	0	0	17	0
SS02CLR4884	2	0	4	6	2	7	0	1	0	92	0	1	3	5	1	0	0	1	1	35	157	1	0	1	1	0	0	1	0	1	0	2	7	0	0	2
SS02CLR7284	468	538	687	932	1,162	1,067	1,131	1,201	1,183	1,082	793	505	498	382	401	752	955	878	999	1,053	1,007	1,283	753	396	397	381	422	742	1,028	1,168	1,216	1,179	1,329	1,161	575	622
SS02LOE7284	0	0	0	0	64	186	176	185	186	191	141	72	92	62	64	139	208	162	188	198	194	217	124	89	71	84	53	140	203	214	221	226	262	239	115	100
SS03ACC7284	30	33	44	48	117	93	133	126	140	101	48	37	29	24	23	40	79	39	57	28	18	15	9	6	17	5	6	8	15	20	21	9	8	19	14	20
SS03BRN6084	39	35	4	0	49	27	45	81	38	49	94	17	26	20	19	37	40	41	6	0	0	0	17	54	19	21	10	30	34	138	59	61	4	0	18	28

Fifteen<- read.csv("OneFive.csv", header = TRUE, stringsAsFactors = FALSE)
Sixteen <- read.csv("OneSix.csv", header = TRUE, stringsAsFactors = FALSE)
FiveSix <- inner_join(Fifteen, Sixteen, by = "Item_Code")
FiveYears <- inner_join(TTF, FiveSix, by = "Item_Code")
kable(head(FiveYears))

Item_Code	X1.1.2012	X2.1.2012	X3.1.2012	X4.1.2012	X5.1.2012	X6.1.2012	X7.1.2012	X8.1.2012	X9.1.2012	X10.1.2012	X11.1.2012	X12.1.2012	X1.1.2013	X2.1.2013	X3.1.2013	X4.1.2013	X5.1.2013	X6.1.2013	X7.1.2013	X8.1.2013	X9.1.2013	X10.1.2013	X11.1.2013	X12.1.2013	X1.1.2014	X2.1.2014	X3.1.2014	X4.1.2014	X5.1.2014	X6.1.2014	X7.1.2014	X8.1.2014	X9.1.2014	X10.1.2014	X11.1.2014	X12.1.2014	Description	X1.1.2015	X2.1.2015	X3.1.2015	X4.1.2015	X5.1.2015	X6.1.2015	X7.1.2015	X8.1.2015	X9.1.2015	X10.1.2015	X11.1.2015	X12.1.2015	X1.12.2016	X2.12.2016	X3.12.2016	X4.12.2016	X5.12.2016	X6.12.2016	X7.12.2016	X8.12.2016	X9.12.2016	X10.12.2016
SS02CLR3660	18	3	5	1	1	2	1	6	0	2	5	6	12	8	3	7	3	2	0	108	1	3	0	3	2	3	1	0	1	0	0	0	0	0	17	0	3/32 CLEAR 36X60	0	0	0	0	0	0	0	0	0	0	1	1	0	0	0	1	0	0	0	0	1	0
SS02CLR4884	2	0	4	6	2	7	0	1	0	92	0	1	3	5	1	0	0	1	1	35	157	1	0	1	1	0	0	1	0	1	0	2	7	0	0	2	3/32 CLEAR 48X84	0	0	0	0	213	0	2	1	0	0	1	0	1	4	7	4	152	2	0	0	1	1
SS02CLR7284	468	538	687	932	1,162	1,067	1,131	1,201	1,183	1,082	793	505	498	382	401	752	955	878	999	1,053	1,007	1,283	753	396	397	381	422	742	1,028	1,168	1,216	1,179	1,329	1,161	575	622	3/32 CLEAR 72 X 84	347	296	566	820	919	1,379	1,159	1,132	1,119	1,238	904	699	512	471	744	789	943	1,341	973	1,255	1,146	1,045
SS02LOE7284	0	0	0	0	64	186	176	185	186	191	141	72	92	62	64	139	208	162	188	198	194	217	124	89	71	84	53	140	203	214	221	226	262	239	115	100	3/32 ENEG ADV LOW-E 72X84	64	45	97	195	192	257	235	250	275	256	208	134	102	82	173	146	232	263	215	262	261	248
SS03ACC7284	30	33	44	48	117	93	133	126	140	101	48	37	29	24	23	40	79	39	57	28	18	15	9	6	17	5	6	8	15	20	21	9	8	19	14	20	1/8 ACCLIMATE 72X84	19	8	27	15	44	26	20	18	20	15	14	10	4	6	8	3	11	15	8	24	9	7
SS03BRN7284	0	0	26	40	0	0	0	0	0	0	0	0	0	0	0	0	0	0	23	42	28	24	6	0	0	0	0	0	0	0	0	0	30	23	7	0	1/8 BRONZE 72X84	0	70	226	23	20	37	32	65	42	36	46	62	59	13	35	59	32	0	0	0	0	0

str(FiveYears)

## 'data.frame':    107 obs. of  60 variables:
##  $ Item_Code  : chr  "SS02CLR3660" "SS02CLR4884" "SS02CLR7284" "SS02LOE7284" ...
##  $ X1.1.2012  : chr  "18" "2" "468" "0" ...
##  $ X2.1.2012  : chr  "3" "0" "538" "0" ...
##  $ X3.1.2012  : chr  "5" "4" "687" "0" ...
##  $ X4.1.2012  : chr  "1" "6" "932" "0" ...
##  $ X5.1.2012  : chr  "1" "2" "1,162" "64" ...
##  $ X6.1.2012  : chr  "2" "7" "1,067" "186" ...
##  $ X7.1.2012  : chr  "1" "0" "1,131" "176" ...
##  $ X8.1.2012  : chr  "6" "1" "1,201" "185" ...
##  $ X9.1.2012  : chr  "0" "0" "1,183" "186" ...
##  $ X10.1.2012 : chr  "2" "92" "1,082" "191" ...
##  $ X11.1.2012 : chr  "5" "0" "793" "141" ...
##  $ X12.1.2012 : chr  "6" "1" "505" "72" ...
##  $ X1.1.2013  : chr  "12" "3" "498" "92" ...
##  $ X2.1.2013  : chr  "8" "5" "382" "62" ...
##  $ X3.1.2013  : chr  "3" "1" "401" "64" ...
##  $ X4.1.2013  : chr  "7" "0" "752" "139" ...
##  $ X5.1.2013  : chr  "3" "0" "955" "208" ...
##  $ X6.1.2013  : chr  "2" "1" "878" "162" ...
##  $ X7.1.2013  : chr  "0" "1" "999" "188" ...
##  $ X8.1.2013  : chr  "108" "35" "1,053" "198" ...
##  $ X9.1.2013  : chr  "1" "157" "1,007" "194" ...
##  $ X10.1.2013 : chr  "3" "1" "1,283" "217" ...
##  $ X11.1.2013 : chr  "0" "0" "753" "124" ...
##  $ X12.1.2013 : chr  "3" "1" "396" "89" ...
##  $ X1.1.2014  : chr  "2" "1" "397" "71" ...
##  $ X2.1.2014  : chr  "3" "0" "381" "84" ...
##  $ X3.1.2014  : chr  "1" "0" "422" "53" ...
##  $ X4.1.2014  : chr  "0" "1" "742" "140" ...
##  $ X5.1.2014  : chr  "1" "0" "1,028" "203" ...
##  $ X6.1.2014  : chr  "0" "1" "1,168" "214" ...
##  $ X7.1.2014  : chr  "0" "0" "1,216" "221" ...
##  $ X8.1.2014  : chr  "0" "2" "1,179" "226" ...
##  $ X9.1.2014  : chr  "0" "7" "1,329" "262" ...
##  $ X10.1.2014 : chr  "0" "0" "1,161" "239" ...
##  $ X11.1.2014 : chr  "17" "0" "575" "115" ...
##  $ X12.1.2014 : chr  "0" "2" "622" "100" ...
##  $ Description: chr  "3/32 CLEAR 36X60" "3/32 CLEAR 48X84" "3/32 CLEAR 72 X 84" "3/32  ENEG ADV LOW-E 72X84" ...
##  $ X1.1.2015  : chr  "0" "0" "347" "64" ...
##  $ X2.1.2015  : chr  "0" "0" "296" "45" ...
##  $ X3.1.2015  : chr  "0" "0" "566" "97" ...
##  $ X4.1.2015  : chr  "0" "0" "820" "195" ...
##  $ X5.1.2015  : chr  "0" "213" "919" "192" ...
##  $ X6.1.2015  : chr  "0" "0" "1,379" "257" ...
##  $ X7.1.2015  : chr  "0" "2" "1,159" "235" ...
##  $ X8.1.2015  : chr  "0" "1" "1,132" "250" ...
##  $ X9.1.2015  : chr  "0" "0" "1,119" "275" ...
##  $ X10.1.2015 : chr  "0" "0" "1,238" "256" ...
##  $ X11.1.2015 : chr  "1" "1" "904" "208" ...
##  $ X12.1.2015 : chr  "1" "0" "699" "134" ...
##  $ X1.12.2016 : chr  "0" "1" "512" "102" ...
##  $ X2.12.2016 : chr  "0" "4" "471" "82" ...
##  $ X3.12.2016 : chr  "0" "7" "744" "173" ...
##  $ X4.12.2016 : chr  "1" "4" "789" "146" ...
##  $ X5.12.2016 : chr  "0" "152" "943" "232" ...
##  $ X6.12.2016 : chr  "0" "2" "1,341" "263" ...
##  $ X7.12.2016 : chr  "0" "0" "973" "215" ...
##  $ X8.12.2016 : chr  "0" "0" "1,255" "262" ...
##  $ X9.12.2016 : chr  "1" "1" "1,146" "261" ...
##  $ X10.12.2016: chr  "0" "1" "1,045" "248" ...

R <- ts(c(34,171,64,48,55,39,33,63,18,45,57,15, 51,33,25,51,34,30,39,34,28,81,78,26, 63,32,36,77,41,31,80,55,72,167,118,142,56,33,187,76,67,73,68,87,113,234,102,108,42,49,99,56,57,42,4,0,112,190), start = c(2012, 1), frequency = 12)
plot(R)

ggseasonplot(R, year.labels=TRUE, continuous=TRUE)

ets(R)

## ETS(A,N,N) 
## 
## Call:
##  ets(y = R) 
## 
##   Smoothing parameters:
##     alpha = 0.1444 
## 
##   Initial states:
##     l = 60.1897 
## 
##   sigma:  46.5554
## 
##      AIC     AICc      BIC 
## 687.0203 687.4647 693.2016

S <- ts(c(91,89,84,77,65,59,73,99,108,167,134,158,191,113,114, 79,46,47,68,78,87,122,90,105,99,109,89,57,69,73,73,60,89,104,109,88,89), start = c(2013, 10), frequency = 12)
plot(S)

ggseasonplot(S, year.labels=TRUE, continuous=TRUE, col = rainbow(12))

ets(S)

## ETS(M,N,N) 
## 
## Call:
##  ets(y = S) 
## 
##   Smoothing parameters:
##     alpha = 0.9999 
## 
##   Initial states:
##     l = 86.3414 
## 
##   sigma:  0.2403
## 
##      AIC     AICc      BIC 
## 365.9519 366.6792 370.7847

SS06BRN96130

U <- ts(c(75,31,42,114,33,53,109,152,71,105,88,51,41,59,29,39,187,79,146,81,59,147,65,71, 48,64,36,46,73,64,94,107,107,108,69,114
, 61,89,48,72,60,52,88,99,114,65,61,95, 71,33,90,41,54,61,83,92,99,120), start = c(2012,1), frequency = 12)
plot(U)

ggseasonplot(U, year.labels=TRUE, continuous=TRUE, col = rainbow(12))

ets(U)

## ETS(A,N,N) 
## 
## Call:
##  ets(y = U) 
## 
##   Smoothing parameters:
##     alpha = 1e-04 
## 
##   Initial states:
##     l = 77.671 
## 
##   sigma:  32.994
## 
##      AIC     AICc      BIC 
## 647.0795 647.5240 653.2609

V <- ts(c(19,26,24,13,18,29,21,37,13,49,30,17, 21,18,22,30,25,30,25,30,24,39,20,27, 13,24,22,29,39,33,35,37,46,49,53,37, 28,37,42,36,40,58,84,77,56,78,66,72,81,77,97,135,82,95,75,85,66,141), start = c(2012,1), frequency = 12)
plot(V)

ggseasonplot(V, year.labels=TRUE, continuous=TRUE, col = rainbow(12), main= "Demand for 3/16 SOLARBAN60 TC")

ets(V)

## ETS(M,A,N) 
## 
## Call:
##  ets(y = V) 
## 
##   Smoothing parameters:
##     alpha = 0.2774 
##     beta  = 1e-04 
## 
##   Initial states:
##     l = 18.1913 
##     b = 1.0596 
## 
##   sigma:  0.3097
## 
##      AIC     AICc      BIC 
## 535.6375 536.7913 545.9397

W <- ts(c(71,56,40, 35,28,36,68,53,74,73,68,85,95,55,45, 30,12,41,42,71,69,79,76,89,110,66,68,31,40,45,57,73,76,92,125,72,68
), start = c(2013, 10), frequency = 12)
plot(W)

ggseasonplot(W,  year.labels=TRUE, continuous=TRUE, col = rainbow(12), main= "Demand for 5/32 ENEG ADV LOW-E 96X130")

ets(W)

## ETS(A,N,N) 
## 
## Call:
##  ets(y = W) 
## 
##   Smoothing parameters:
##     alpha = 0.8778 
## 
##   Initial states:
##     l = 68.9324 
## 
##   sigma:  20.2136
## 
##      AIC     AICc      BIC 
## 362.0743 362.8016 366.9071

SS10STP96130

X <- ts(c(17,33,24,25,17,20,14,14,12,11,15,10
,16,13,18,26,41,16,5,33,47,36,20,24
, 19,25,56,30,48,47,40,24,31,41,35,40, 28,29,32,49,63,84,55,70,70,92,62,49,42,69,79,71,85,72,93,52,43), start = c(2012,1), frequency = 12)
plot(X)

ggseasonplot(X, year.labels=TRUE, continuous=TRUE, col = rainbow(12), main= "Demand for 3/8 STARPHIRE 96X130")

ets(X)

## ETS(A,N,N) 
## 
## Call:
##  ets(y = X) 
## 
##   Smoothing parameters:
##     alpha = 0.4705 
## 
##   Initial states:
##     l = 22.4425 
## 
##   sigma:  14.1105
## 
##      AIC     AICc      BIC 
## 538.2029 538.6557 544.3320

SS12STP96130

Y <- ts(c(12,11,17,17,43,20,17,23,12,18,11,29, 13,30,17,6,21,35,28,26,18,101,52,18, 23,23,19,39,35,27,47,39,54,38,34,34, 
12,32,27,43,56,60,51,31,40,52,39,29,47,42,62,54,107,81,96,112,110,60), start = c(2012,12), frequency = 12)
plot(Y)

ggseasonplot(Y, year.labels=TRUE, continuous=TRUE, col = rainbow(12), main= "1/2 STARPHIRE 96X130")

ets(Y)

## ETS(M,N,A) 
## 
## Call:
##  ets(y = Y) 
## 
##   Smoothing parameters:
##     alpha = 0.2394 
##     gamma = 1e-04 
## 
##   Initial states:
##     l = 21.8664 
##     s=0.6057 -4.8189 23.0302 -7.6226 -3.226 3.0836
##            1.1993 8.0777 -5.0708 -5.39 -2.0166 -7.8515
## 
##   sigma:  0.4419
## 
##      AIC     AICc      BIC 
## 567.0208 578.4494 597.9275

Z <- ts(c(26,88,31,85,33,28,81,60,38,26,35,32, 15,17,65,84,102,33,12,48,133,402,34,53, 90,17,25,22,135,52,76,11,39,72,155,78, 30,5,93,245,34,56,56,45,40,2,0,0,28,19,92,54,55,98,63,16,48,63), start = c(2012, 1), frequency = 12)
plot(Z)

ggseasonplot(Z, year.labels=TRUE, continuous=TRUE, col = rainbow(12), main= "1/4 SN-68 TC 96X130")

ets(Z)

## ETS(A,N,N) 
## 
## Call:
##  ets(y = Z) 
## 
##   Smoothing parameters:
##     alpha = 1e-04 
## 
##   Initial states:
##     l = 59.9237 
## 
##   sigma:  62.0443
## 
##      AIC     AICc      BIC 
## 720.3362 720.7806 726.5175

AA <- ts(c(12,16,34,17,13,7,16,7,46,48,34,16
,25,50,14,32,39,26,47,61,30,44,26,7
,5,30,18,28,49,44,48,70,87,66,45,55
, 53,30,29,44,68,63,67,56,73,54,28,48,39,64,36,22,55,91,44,44,10,24), start = c(2012,1), frequency = 12)
plot(AA)

ggseasonplot(AA, year.labels=TRUE, continuous=TRUE, col = rainbow(12), main= "1/4 CLEAR LAMI 40X80")

ets(AA)

## ETS(A,N,N) 
## 
## Call:
##  ets(y = AA) 
## 
##   Smoothing parameters:
##     alpha = 0.4678 
## 
##   Initial states:
##     l = 16.4249 
## 
##   sigma:  18.2727
## 
##      AIC     AICc      BIC 
## 578.5333 578.9777 584.7146

BB <- ts(c(55,19,11,25,15,51,9,9,44,65,21,22
, 24,11,18,4,35,22,26,29,26,30,8,24
, 29,7,14,9,12,12,9,15,7,10,22,14
, 39,16,16,42,80,38,24,24,31,24,40,166,50,30,47,74,24,11,25,27,46,63), start = c(2012,1), frequency = 12)
plot(BB)

ggseasonplot(BB, year.labels=TRUE, continuous=TRUE, col = rainbow(12), main= "1/4 GREEN 96X130
")

ets(BB)

## ETS(M,N,M) 
## 
## Call:
##  ets(y = BB) 
## 
##   Smoothing parameters:
##     alpha = 0.275 
##     gamma = 1e-04 
## 
##   Initial states:
##     l = 26.25 
##     s=1.7786 0.8197 1.1797 0.9376 0.6391 0.6134
##            1.0326 1.515 0.9358 0.581 0.5011 1.4664
## 
##   sigma:  0.5803
## 
##      AIC     AICc      BIC 
## 575.1591 586.5877 606.0658

CC <- ts(c(12,31,10,14,25,25,12,34,19,26,9,13
, 12,44,25,17,19,26,37,31,13,27,8,23
, 11,8,9,18,32,36,30,42,24,14,20,1,
5,9,23,8,49,54,58,47,49,103,26,47,33,35,32,35,45,44,67,74,44,59), start = c(2012,1), frequency = 12)
plot(CC)

ggseasonplot(CC, year.labels=TRUE, continuous=TRUE, col = rainbow(12), main= "3/16 SOLARBAN60 96X130")

ets(CC)

## ETS(M,N,A) 
## 
## Call:
##  ets(y = CC) 
## 
##   Smoothing parameters:
##     alpha = 0.1806 
##     gamma = 1e-04 
## 
##   Initial states:
##     l = 20.6524 
##     s=-5.9888 -9.5895 12.3898 -3.7328 7.6394 4.7185
##            6.7458 7.9288 -7.3201 -3.9303 2.4739 -11.3347
## 
##   sigma:  0.4641
## 
##      AIC     AICc      BIC 
## 543.1273 554.5559 574.0340

DD <- ts(c(15,25,22,20,30,35,30,48,27,23,25,13
,23,17,26,35,34,39,41,29,27,25,29,19
,20,20,21,35,32,52,58,42,42,30,30,39
,25,26,48,47,33,56,48,58,81,48,56,30,15,24,30,16,16,39,41,58,53,46), start = c(2012,1), frequency = 12)
plot(DD)

ggseasonplot(DD, year.labels=TRUE, continuous=TRUE, col = rainbow(12), main= "1/2 SHOWERGUARD 96X130")

ets(DD)

## ETS(M,N,M) 
## 
## Call:
##  ets(y = DD) 
## 
##   Smoothing parameters:
##     alpha = 0.171 
##     gamma = 1e-04 
## 
##   Initial states:
##     l = 28.0393 
##     s=0.7311 0.9579 0.9549 1.3279 1.3438 1.2214
##            1.3056 0.9407 0.9978 0.8951 0.6837 0.64
## 
##   sigma:  0.263
## 
##      AIC     AICc      BIC 
## 510.2532 521.6818 541.1599

EE<- ts(c(32,16,21,23,41,72,42,37,62,32,40,33
, 54,28,57,28,31,44,25,27,38,32,31,33
, 32,19,24,21,31,27,44,40,50,38,21,32, 17,19,28,21,31,25,56,32,49,36,34,23,25,34,39,39,74,65,40,47,35
), start = c(2012,1), frequency = 12)
plot(EE)

ggseasonplot(EE, year.labels=TRUE, continuous=TRUE, col = rainbow(12), main= "3/16 BRONZE 96X130")

ets(EE)

## ETS(M,N,N) 
## 
## Call:
##  ets(y = EE) 
## 
##   Smoothing parameters:
##     alpha = 1e-04 
## 
##   Initial states:
##     l = 35.5607 
## 
##   sigma:  0.3692
## 
##      AIC     AICc      BIC 
## 529.9763 530.4291 536.1054

FF <- ts(c(14,6,13,23,18,14,11,47,72,21,23,17
,19,16,19,14,32,18,19,17,23,34,36,12
,24,14,26,22,18,19,35,24,37,57,33,68
,34,24,35,16,17,58,28,38,20,56,17,22,16,18,22,23,18,78,56,117,36,43), start = c(2012,1), frequency = 12)
plot(FF)

ggseasonplot(FF, year.labels=TRUE, continuous=TRUE, col = rainbow(12), main= "1/4 ACID ETCH 96X130")

ets(FF)

## ETS(M,N,M) 
## 
## Call:
##  ets(y = FF) 
## 
##   Smoothing parameters:
##     alpha = 0.1061 
##     gamma = 1e-04 
## 
##   Initial states:
##     l = 22.5386 
##     s=1.0281 0.9706 1.2881 1.7521 1.5068 0.9932
##            1.1642 0.7424 0.6114 0.807 0.4764 0.6597
## 
##   sigma:  0.5026
## 
##      AIC     AICc      BIC 
## 559.9694 571.3980 590.8761

GG <- ts(c(17,23,32,65,16,33,8,11,40,36,27,25
,16,39,22,72,25,14,32,39,32,20,9,32
,9,27,16,12,31,18,41,21,30,56,25,37
,8,18,28,34,23,20,44,36,64,44,28,28,38,40,36,29,43,47,27,45,24,15), start = c(2012,1), frequency = 12)
plot(GG)

ggseasonplot(GG, year.labels=TRUE, continuous=TRUE, col = rainbow(12), main= "1/4 CLEAR LAMI 40X80")

ets(GG)

## ETS(A,N,N) 
## 
## Call:
##  ets(y = GG) 
## 
##   Smoothing parameters:
##     alpha = 1e-04 
## 
##   Initial states:
##     l = 29.7751 
## 
##   sigma:  13.9551
## 
##      AIC     AICc      BIC 
## 547.2638 547.7082 553.4451

1/4 SOLARBAN 70XL 100X144

AA <- ts(c(18,35,9,50,19,7,10,6,14,8,8,36,3,17,8,62,61,10,148,17,95
, 12,14,20,6,22,16,20,17,16,29,12,28,19,41,24,37,33,20,6,59,23,40), start = c(2012,4), frequency = 12)
plot(AA)

ggseasonplot(AA, year.labels=TRUE, continuous=TRUE, col = rainbow(12), main= "1/4 SOLARBAN 70XL 100X144")

ets(AA)

## ETS(M,N,A) 
## 
## Call:
##  ets(y = AA) 
## 
##   Smoothing parameters:
##     alpha = 0.041 
##     gamma = 1e-04 
## 
##   Initial states:
##     l = 24.9302 
##     s=-0.6091 -2.3216 -14.0505 17.0458 -16.6648 40.4515
##            -13.1127 12.2626 9.1794 -15.0913 -7.3387 -9.7505
## 
##   sigma:  0.6198
## 
##      AIC     AICc      BIC 
## 413.4897 431.2675 439.9077

=C2&“,”&D2&“,”&E2&“,”&F2&“,”&G2&“,”&H2&“,”&I2&“,”&J2&“,”&K2&“,”&L2&“,”&M2&“,”&N2&“,”&O2&“,”&P2&“,”&Q2&“,”&R2&“,”&S2&“,”&T2&“,”&U2&“,”&V2&“,”&W2&“,”&X2&“,”&Y2&“,”&Z2&“,”&AA2&“,”&AB2&“,”&AC2&“,”&AD2&“,”&AE2&“,”&AF2&“,”&AG2&“,”&AH2&“,”&AI2&“,”&AJ2&“,”&AK2&“,”&AL2&“,”&AM2&“,”&AN2&“,”&AO2&“,”&AP2&“,”&AQ2&“,”&AR2&“,”&AS2&“,”&AT2&“,”&AU2&“,”&AV2&“,”&AW2&“,”&AX2&“,”&AY2&“,”&AZ2&“,”&BA2&“,”&BB2&“,”&BC2&“,”&BD2&“,”&BE2&“,”&BF2&“,”&BG2&“,”&BH2

Group <- cbind(A,B,C,D,E,FF,G,H,I,J,K,L,M,N,O,P,Q) kable(head(Group)) RowGroup <- rbind(A,B,C,D,E,FF,G,H,I,J,K,L,M,N,O,P,Q) kable(head(RowGroup))