Pakistan Solar Panel Prices and Global Polysilicon: An ARIMAX Analysis

Cointegration, Granger causality, and a 12-month forecast, with policy notes

Author

Prof. Dr. Zahid Asghar

Published

September 22, 2026

1 Overview

This note models monthly Pakistan solar panel prices (log scale) against global polysilicon prices (log scale), January 2017 to May 2026 (n = 113). Because panel prices in Pakistan are largely a function of an imported input, the analysis treats polysilicon as an exogenous driver rather than modelling panel price in isolation, and checks whether the two series are cointegrated before forecasting.

Code

pkgs <- c("forecast", "tseries", "urca", "lmtest", "ggplot2")
invisible(lapply(pkgs, function(p) if (!requireNamespace(p, quietly = TRUE))
  install.packages(p, repos = "https://cran.r-project.org")))
library(forecast); library(tseries); library(urca); library(lmtest); library(ggplot2)

brand_green <- "#2C5F2D"
brand_brown <- "#8B7355"

2 Data preparation

Code

df <- read.csv(params$data_path, stringsAsFactors = FALSE)
df$Date <- as.Date(df$Date, format = "%d/%m/%Y")
df <- df[order(df$Date), ]
stopifnot(!anyNA(df$log_commodity_price), !anyNA(df$log_polysilicon))

y <- ts(df$log_commodity_price, start = c(2017, 1), frequency = 12)  # panel price (log)
x <- ts(df$log_polysilicon,     start = c(2017, 1), frequency = 12)  # polysilicon price (log)

knitr::kable(head(df, 6), caption = "First six observations")

First six observations
Area	Date	Commodity.type	log_commodity_price	log_polysilicon
Pakistan	2017-01-01	Panels	-0.97	2.86
Pakistan	2017-02-01	Panels	-0.97	2.87
Pakistan	2017-03-01	Panels	-0.97	2.81
Pakistan	2017-04-01	Panels	-0.97	2.62
Pakistan	2017-05-01	Panels	-0.97	2.68
Pakistan	2017-06-01	Panels	-0.97	2.71

3 Exploratory analysis

Code

autoplot(cbind(PanelPrice = y, Polysilicon = x), facets = TRUE) +
  scale_color_manual(values = c(brand_green, brand_brown)) +
  labs(x = NULL, y = NULL) + theme_minimal()

Log panel price (Pakistan) vs log polysilicon price

Code

stl_y <- stl(y, s.window = "periodic")
autoplot(stl_y) + theme_minimal()

4 Stationarity tests

Both series look non-stationary in levels. The battery below runs ADF, KPSS, and Phillips-Perron on levels and first differences.

Code

unit_root_battery <- function(series, name) {
  cat("\n===", name, "===\n")
  cat("ADF (level):  p =", round(adf.test(series)$p.value, 4), "\n")
  cat("KPSS (level): p =", round(kpss.test(series, "Level")$p.value, 4), "\n")
  cat("PP  (level):  p =", round(pp.test(series)$p.value, 4), "\n")
  cat("ADF (1st diff): p =", round(adf.test(diff(series))$p.value, 4), "\n")
}
unit_root_battery(y, "Log Panel Price")


=== Log Panel Price ===
ADF (level):  p = 0.7363

KPSS (level): p = 0.01 
PP  (level):  p = 0.846 
ADF (1st diff): p = 0.0238

Code

unit_root_battery(x, "Log Polysilicon Price")


=== Log Polysilicon Price ===
ADF (level):  p = 0.8965 
KPSS (level): p = 0.0102 
PP  (level):  p = 0.9251 
ADF (1st diff): p = 0.1541

Code

cat("\nndiffs(y) =", ndiffs(y), " | ndiffs(x) =", ndiffs(x), "\n")


ndiffs(y) = 1  | ndiffs(x) = 1

5 ACF / PACF identification

Code

par(mfrow = c(2, 2))
acf(diff(y),  main = "ACF: d(Panel Price)")
pacf(diff(y), main = "PACF: d(Panel Price)")
acf(diff(x),  main = "ACF: d(Polysilicon)")
pacf(diff(x), main = "PACF: d(Polysilicon)")

Code

par(mfrow = c(1, 1))

6 Cointegration

Code

eg_reg <- lm(y ~ x)
summary(eg_reg)


Call:
lm(formula = y ~ x)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.47923 -0.23939 -0.01016  0.21182  0.52348 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -2.75081    0.08884  -30.96   <2e-16 ***
x            0.47990    0.03493   13.74   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.263 on 111 degrees of freedom
Multiple R-squared:  0.6297,    Adjusted R-squared:  0.6263 
F-statistic: 188.7 on 1 and 111 DF,  p-value: < 2.2e-16

Code

eg_resid <- residuals(eg_reg)
eg_adf <- ur.df(eg_resid, type = "none", selectlags = "AIC")
cat("\nEngle-Granger step 2 (residual ADF), tau =", round(eg_adf@teststat[1], 3),
    " | 5% crit =", eg_adf@cval[1, "5pct"], "\n")


Engle-Granger step 2 (residual ADF), tau = -1.975  | 5% crit = -1.95

Code

jo <- ca.jo(cbind(y, x), type = "trace", ecdet = "const", K = 2)
summary(jo)


###################### 
# Johansen-Procedure # 
###################### 

Test type: trace statistic , without linear trend and constant in cointegration 

Eigenvalues (lambda):
[1] 4.429505e-02 4.028875e-02 4.163336e-17

Values of teststatistic and critical values of test:

         test 10pct  5pct  1pct
r <= 1 | 4.56  7.52  9.24 12.97
r = 0  | 9.59 17.85 19.96 24.60

Eigenvectors, normalised to first column:
(These are the cointegration relations)

               y.l2       x.l2   constant
y.l2      1.0000000  1.0000000 1.00000000
x.l2     -0.5163466 -0.9587116 0.47322825
constant  2.9994605  3.4414821 0.07950448

Weights W:
(This is the loading matrix)

           y.l2       x.l2      constant
y.d -0.02888688 0.01379029 -2.245622e-17
x.d  0.05266642 0.02838822  3.825622e-17

Sample-size caveat

With n = 113 monthly observations, the Johansen and Engle-Granger tests have limited power. Read the cointegration result as indicative and re-estimate as more months of data accrue.

7 Granger causality

Code

cat("Does polysilicon Granger-cause panel price?\n")

Does polysilicon Granger-cause panel price?

Code

print(grangertest(y ~ x, order = 3))

Granger causality test

Model 1: y ~ Lags(y, 1:3) + Lags(x, 1:3)
Model 2: y ~ Lags(y, 1:3)
  Res.Df Df      F Pr(>F)
1    103                 
2    106 -3 0.5479 0.6507

Code

cat("\nDoes panel price Granger-cause polysilicon?\n")


Does panel price Granger-cause polysilicon?

Code

print(grangertest(x ~ y, order = 3))

Granger causality test

Model 1: x ~ Lags(x, 1:3) + Lags(y, 1:3)
Model 2: x ~ Lags(x, 1:3)
  Res.Df Df     F   Pr(>F)   
1    103                     
2    106 -3 4.583 0.004708 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

8 Model estimation and out-of-sample accuracy

The last 13 months are held out; ARIMA, ARIMAX, and a naive random-walk benchmark are compared on RMSE, MAE, and MAPE.

Code

train_h  <- 100
y_train  <- window(y, end = time(y)[train_h])
x_train  <- window(x, end = time(x)[train_h])
y_test   <- window(y, start = time(y)[train_h + 1])
x_test   <- window(x, start = time(x)[train_h + 1])

m_arima  <- auto.arima(y_train, stepwise = FALSE, approximation = FALSE)
m_arimax <- auto.arima(y_train, xreg = x_train, stepwise = FALSE, approximation = FALSE)

summary(m_arima)

Series: y_train 
ARIMA(0,1,0)(0,0,2)[12] with drift 

Coefficients:
       sma1    sma2    drift
      0.142  0.2932  -0.0145
s.e.  0.104  0.1249   0.0076

sigma^2 = 0.003246:  log likelihood = 143.55
AIC=-279.1   AICc=-278.68   BIC=-268.72

Training set error measures:
                        ME       RMSE        MAE        MPE     MAPE      MASE
Training set -2.328004e-05 0.05582143 0.03535319 -0.1112613 2.349967 0.1578428
                    ACF1
Training set -0.08212039

Code

summary(m_arimax)

Series: y_train 
Regression with ARIMA(1,1,0)(0,0,2)[12] errors 

Coefficients:
          ar1    sma1    sma2    drift    xreg
      -0.1588  0.1308  0.3371  -0.0131  0.0875
s.e.   0.1022  0.1032  0.1323   0.0066  0.0386

sigma^2 = 0.00312:  log likelihood = 146.19
AIC=-280.37   AICc=-279.46   BIC=-264.8

Training set error measures:
                        ME       RMSE        MAE        MPE     MAPE      MASE
Training set -5.372398e-06 0.05415472 0.03418746 -0.1111063 2.290322 0.1526381
                    ACF1
Training set 0.007790024

Code

checkresiduals(m_arimax)


    Ljung-Box test

data:  Residuals from Regression with ARIMA(1,1,0)(0,0,2)[12] errors
Q* = 19.171, df = 17, p-value = 0.3188

Model df: 3.   Total lags used: 20

Code

cat("Jarque-Bera (ARIMAX resid): p =",
    round(jarque.bera.test(residuals(m_arimax))$p.value, 4), "\n")

Jarque-Bera (ARIMAX resid): p = 0

Code

cat("ARCH-LM (ARIMAX resid, lag 6):\n")

ARCH-LM (ARIMAX resid, lag 6):

Code

print(Box.test(residuals(m_arimax)^2, lag = 6, type = "Ljung-Box"))


    Box-Ljung test

data:  residuals(m_arimax)^2
X-squared = 2.0366, df = 6, p-value = 0.9163

Code

h <- length(y_test)
f_arima  <- forecast(m_arima,  h = h)
f_arimax <- forecast(m_arimax, xreg = x_test, h = h)   # oracle: actual future polysilicon
f_naive  <- naive(y_train, h = h)

acc_tbl <- rbind(
  ARIMA  = accuracy(f_arima,  y_test)["Test set", c("RMSE", "MAE", "MAPE")],
  ARIMAX = accuracy(f_arimax, y_test)["Test set", c("RMSE", "MAE", "MAPE")],
  Naive  = accuracy(f_naive,  y_test)["Test set", c("RMSE", "MAE", "MAPE")]
)
knitr::kable(round(acc_tbl, 4), caption = "Out-of-sample accuracy, 13-month holdout")

Out-of-sample accuracy, 13-month holdout
	RMSE	MAE	MAPE
ARIMA	0.2819	0.2388	10.7452
ARIMAX	0.2487	0.2052	9.2516
Naive	0.0900	0.0585	2.6350

Code

cv_arima <- tsCV(y, function(z, h) forecast(auto.arima(z), h = h), h = 1)
rmse_cv  <- round(sqrt(mean(cv_arima^2, na.rm = TRUE)), 4)
cat("Rolling 1-step CV RMSE (ARIMA):", rmse_cv, "\n")

Rolling 1-step CV RMSE (ARIMA): 0.0642

9 Final model and 12-month forecast

Code

final_arimax <- auto.arima(y, xreg = x, stepwise = FALSE, approximation = FALSE)
final_x_ar   <- auto.arima(x)
x_future     <- forecast(final_x_ar, h = 12)$mean

fc_trend <- forecast(final_arimax, xreg = x_future, h = 12)
fc_flat  <- forecast(final_arimax, xreg = rep(tail(x, 1), 12), h = 12)

autoplot(fc_trend) +
  labs(subtitle = "xreg = ARIMA-projected polysilicon price", y = "log(panel price)") +
  theme_minimal()

12-month forecast: log panel price (Pakistan)

Code

fc_table <- data.frame(
  Month          = seq(as.Date("2026-06-01"), by = "month", length.out = 12),
  Trend_scenario = as.numeric(fc_trend$mean),
  Flat_scenario  = as.numeric(fc_flat$mean),
  Lo95           = as.numeric(fc_trend$lower[, 2]),
  Hi95           = as.numeric(fc_trend$upper[, 2])
)
knitr::kable(round(fc_table[, -1], 3), caption = "12-month forecast (log scale)")

12-month forecast (log scale)
Trend_scenario	Flat_scenario	Lo95	Hi95
-2.104	-2.106	-2.218	-1.990
-2.107	-2.109	-2.254	-1.961
-2.116	-2.118	-2.308	-1.925
-2.125	-2.127	-2.353	-1.897
-2.134	-2.136	-2.394	-1.875
-2.143	-2.145	-2.431	-1.856
-2.152	-2.154	-2.465	-1.839
-2.161	-2.163	-2.498	-1.825
-2.170	-2.172	-2.529	-1.811
-2.179	-2.181	-2.558	-1.800
-2.188	-2.190	-2.587	-1.789
-2.197	-2.199	-2.615	-1.779

Code

trend_dir <- ifelse(tail(fc_table$Trend_scenario, 1) < fc_table$Trend_scenario[1],
                     "continued decline", "stabilization")

10 Policy relevance

Policy relevance

Pass-through: the cointegrating regression implies a 1% rise in global polysilicon price is associated with roughly a 0.48% long-run change in Pakistan panel prices.
Direction of influence: the Granger results point to the global polysilicon market driving domestic panel prices rather than the reverse, consistent with Pakistan being a price-taker with negligible domestic module manufacturing.
Forecast signal: the 12-month projection points to a continued decline in log panel prices, which is supportive of further residential and commercial solar uptake so long as duty and tax structures on imported modules are not tightened.
Levers with real traction: since pricing is driven by an imported input with no domestic hedge, the policy options with the most bite are
1. tariff/GST rationalization on panels and inverters, (b) exchange-rate pass-through management, and (c) incentives for local cell/module assembly to reduce exposure to polysilicon price cycles.

11 Reproducibility

Code

sessionInfo()

R version 4.6.0 (2026-04-24 ucrt)
Platform: x86_64-w64-mingw32/x64
Running under: Windows 11 x64 (build 22621)

Matrix products: default
  LAPACK version 3.12.1

locale:
[1] LC_COLLATE=English_Pakistan.utf8  LC_CTYPE=English_Pakistan.utf8   
[3] LC_MONETARY=English_Pakistan.utf8 LC_NUMERIC=C                     
[5] LC_TIME=English_Pakistan.utf8    

time zone: Asia/Karachi
tzcode source: internal

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] ggplot2_4.0.3   lmtest_0.9-40   zoo_1.8-15      urca_1.3-4     
[5] tseries_0.10-61 forecast_9.0.2 

loaded via a namespace (and not attached):
 [1] generics_0.1.4     lattice_0.22-9     digest_0.6.39      magrittr_2.0.5    
 [5] evaluate_1.0.5     grid_4.6.0         RColorBrewer_1.1-3 fastmap_1.2.0     
 [9] jsonlite_2.0.0     scales_1.4.0       cli_3.6.6          rlang_1.2.0       
[13] withr_3.0.2        yaml_2.3.12        otel_0.2.0         tools_4.6.0       
[17] parallel_4.6.0     dplyr_1.2.1        colorspace_2.1-2   curl_7.1.0        
[21] vctrs_0.7.3        R6_2.6.1           stats4_4.6.0       lifecycle_1.0.5   
[25] htmlwidgets_1.6.4  pkgconfig_2.0.3    pillar_1.11.1      gtable_0.3.6      
[29] glue_1.8.1         quantmod_0.4.28    Rcpp_1.1.1-1       xfun_0.57         
[33] tibble_3.3.1       tidyselect_1.2.1   rstudioapi_0.18.0  knitr_1.51        
[37] farver_2.1.2       htmltools_0.5.9    nlme_3.1-169       rmarkdown_2.31    
[41] xts_0.14.2         labeling_0.4.3     timeDate_4052.112  fracdiff_1.5-3    
[45] compiler_4.6.0     S7_0.2.2           quadprog_1.5-8     TTR_0.24.4