Piotroski F-score and Future Stock Returns: Evidence from Taiwan’s Textile Industry
Business Valuation Project — Econometric Output
Angelina (蔣舜涵) — Student ID 111034049
2026-06-13
1 Data
Input: panel_clean.csv — merged quarterly panel for
listed Taiwanese textile firms (TSE Textiles & Fiber), assembled
from five TEJ Pro exports (balance sheet, income statement, cash-flow
statement, stock returns, per-share book value). Net income =
consolidated total profit after tax; all flow items are single-quarter
figures.
d <- read.csv("panel_clean.csv", )
d$t <- d$year * 4 + d$quarter
d$yq <- paste0(d$year, "Q", d$quarter)
d <- d[order(d$coid, d$t), ]
cat("Observations:", nrow(d), " Firms:", length(unique(d$coid)),
" Period:", min(d$yq), "->", max(d$yq), "\n")## Observations: 1954 Firms: 42 Period: 2014Q2 -> 2026Q12 Variable Construction
# continuous financial ratios
d$roa <- d$ni / d$ta
d$cfo_ta <- d$cfo / d$ta
d$cr <- d$ca / d$cl
d$gm <- (d$rev - d$cogs) / d$rev
d$ato <- d$rev / d$ta
d$ltd_r <- d$ltd / d$ta
d$debt_r <- d$tl / d$ta
d$size <- log(d$ta) # firm size = log(total assets)
d$pb <- d$price / d$bvps # price-to-book
# year-over-year change (lag 4 quarters) via explicit (coid, t-4) merge
lagvars <- c("roa","cr","gm","ato","ltd_r","cs")
lg <- d[, c("coid","t",lagvars)]; lg$t <- lg$t + 4
names(lg)[-(1:2)] <- paste0(lagvars,"_l4")
d <- merge(d, lg, by = c("coid","t"), all.x = TRUE)
d <- d[order(d$coid, d$t), ]
d$d_roa <- d$roa - d$roa_l4
d$d_cr <- d$cr - d$cr_l4
d$d_gm <- d$gm - d$gm_l4
d$d_ato <- d$ato - d$ato_l4
d$d_ltd <- d$ltd_r - d$ltd_r_l4
# 9 Piotroski F-score signals (NA where inputs unavailable)
d$F1 <- as.integer(d$roa > 0)
d$F2 <- as.integer(d$cfo > 0)
d$F3 <- as.integer(d$d_roa > 0)
d$F4 <- as.integer(d$cfo_ta > d$roa)
d$F5 <- as.integer(d$d_ltd < 0)
d$F6 <- as.integer(d$d_cr > 0)
d$F7 <- as.integer(d$cs <= d$cs_l4)
d$F8 <- as.integer(d$d_gm > 0)
d$F9 <- as.integer(d$d_ato > 0)
d$Fscore <- rowSums(d[, paste0("F",1:9)]) # NA if any signal is NA
# future one-quarter return via (coid, t-1) merge
fu <- d[, c("coid","t","ret")]; fu$t <- fu$t - 1
names(fu)[3] <- "fut_ret"
d <- merge(d, fu, by = c("coid","t"), all.x = TRUE)
d <- d[order(d$coid, d$t), ]
# Winner = future return above industry (same-quarter) median
med <- ave(d$fut_ret, d$yq, FUN = function(x) median(x, na.rm = TRUE))
d$Winner <- as.integer(d$fut_ret > med)
# five performance classes = within-quarter quintiles of future return
d$Cat5 <- ave(d$fut_ret, d$yq, FUN = function(x){
out <- rep(NA_integer_, length(x))
ok <- is.finite(x)
if (sum(ok) < 5) return(out)
b <- unique(quantile(x[ok], seq(0,1,.2), na.rm = TRUE))
if (length(b) < 2) return(out)
out[ok] <- as.integer(cut(x[ok], b, include.lowest = TRUE))
out
})
reg <- c("F1","F2","F3","F4","F5","F6","F7","F8","F9","Fscore","fut_ret",
"roa","cfo_ta","d_roa","d_cr","d_gm","d_ato","size","pb","debt_r","Winner")
D <- d[complete.cases(d[, reg]), ]
cat("Analysis sample N =", nrow(D), " (", min(D$yq), "->", max(D$yq), ")\n")## Analysis sample N = 1691 ( 2015Q2 -> 2025Q3 )3 Descriptive Statistics
dv <- c("Fscore","roa","cfo_ta","d_roa","d_cr","d_gm","d_ato","fut_ret","size","pb","debt_r")
round(t(sapply(D[,dv], function(x) c(N=length(x), Mean=mean(x), SD=sd(x),
Min=min(x), Median=median(x), Max=max(x)))), 4)## N Mean SD Min Median Max
## Fscore 1691 5.1508 1.7733 1.0000 5.0000 9.0000
## roa 1691 0.0060 0.0338 -0.2440 0.0039 0.7209
## cfo_ta 1691 0.0107 0.0424 -0.2190 0.0078 0.4340
## d_roa 1691 -0.0003 0.0449 -0.7202 -0.0006 0.9649
## d_cr 1691 0.0104 5.3780 -62.3026 -0.0016 75.1443
## d_gm 1691 -0.0014 0.2797 -4.9803 0.0007 5.2988
## d_ato 1691 -0.0063 0.0557 -0.7021 -0.0031 0.4334
## fut_ret 1691 1.4253 18.3528 -61.0990 -0.3635 250.0007
## size 1691 15.7747 1.3756 10.9641 15.4549 20.3406
## pb 1691 1.4600 1.3640 0.1843 1.0151 16.7255
## debt_r 1691 0.4150 0.1991 0.0329 0.4156 0.9843cat("\nF-score distribution:\n"); print(table(D$Fscore))##
## F-score distribution:##
## 1 2 3 4 5 6 7 8 9
## 25 101 198 271 368 325 246 120 374 Correlation Analysis
cv <- c("Fscore","roa","cfo_ta","d_roa","d_cr","d_gm","d_ato","fut_ret")
round(cor(D[,cv]), 3)## Fscore roa cfo_ta d_roa d_cr d_gm d_ato fut_ret
## Fscore 1.000 0.272 0.381 0.218 0.009 0.215 0.337 0.053
## roa 0.272 1.000 0.199 0.663 -0.028 0.178 0.117 -0.003
## cfo_ta 0.381 0.199 1.000 0.007 -0.027 0.043 -0.044 0.001
## d_roa 0.218 0.663 0.007 1.000 -0.034 0.189 0.252 0.011
## d_cr 0.009 -0.028 -0.027 -0.034 1.000 -0.223 -0.048 0.013
## d_gm 0.215 0.178 0.043 0.189 -0.223 1.000 0.145 -0.014
## d_ato 0.337 0.117 -0.044 0.252 -0.048 0.145 1.000 -0.005
## fut_ret 0.053 -0.003 0.001 0.011 0.013 -0.014 -0.005 1.0005 OLS Regression
m0 <- lm(fut_ret ~ Fscore, data = D) # simple
m1 <- lm(fut_ret ~ F1+F2+F3+F4+F5+F6+F7+F8+F9, data = D) # Model 1
m2 <- lm(fut_ret ~ roa+cfo_ta+d_roa+d_cr+d_gm+d_ato, data = D) # Model 2
cat("Simple OLS (FutRet ~ F-score):\n"); print(summary(m0)$coefficients)## Simple OLS (FutRet ~ F-score):## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.41708 1.36982 -1.0345 0.301049
## Fscore 0.55184 0.25147 2.1945 0.028336cat(sprintf("R2 = %.4f p = %.4f\n\n", summary(m0)$r.squared, summary(m0)$coefficients[2,4]))## R2 = 0.0028 p = 0.0283cat("Model 1 (F1-F9): adj R2 = %.4f\n", summary(m1)$adj.r.squared)## Model 1 (F1-F9): adj R2 = %.4f
## 0.0028117print(round(summary(m1)$coefficients,4))## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.6698 1.6681 -1.0010 0.3170
## F1 0.4451 1.0785 0.4127 0.6799
## F2 -0.5950 1.2973 -0.4586 0.6466
## F3 2.5851 1.0611 2.4363 0.0149
## F4 1.3082 1.2359 1.0585 0.2900
## F5 0.6639 0.9390 0.7071 0.4796
## F6 -0.3346 0.9350 -0.3579 0.7205
## F7 1.4088 1.2691 1.1101 0.2671
## F8 0.7919 1.0098 0.7842 0.4330
## F9 -1.3255 0.9473 -1.3993 0.1619cat("\nModel 2 (continuous): adj R2 =", round(summary(m2)$adj.r.squared,4), "\n")##
## Model 2 (continuous): adj R2 = -0.0029print(round(summary(m2)$coefficients,4))## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.4542 0.4703 3.0922 0.0020
## roa -10.6642 18.3522 -0.5811 0.5613
## cfo_ta 2.1008 10.9252 0.1923 0.8475
## d_roa 11.5939 13.9050 0.8338 0.4045
## d_cr 0.0346 0.0853 0.4058 0.6849
## d_gm -0.8168 1.6824 -0.4855 0.6274
## d_ato -2.3718 8.3764 -0.2831 0.77716 Logistic Regression
acc <- function(m) mean((predict(m,type="response")>.5) == m$model[[1]])
m3a <- glm(Winner ~ F1+F2+F3+F4+F5+F6+F7+F8+F9, data=D, family=binomial)
m3b <- glm(Winner ~ Fscore+size+pb+roa+debt_r, data=D, family=binomial)
m4 <- glm(Winner ~ roa+cfo_ta+d_roa+d_cr+d_gm+d_ato, data=D, family=binomial)
cat(sprintf("Accuracy 3A=%.3f 3B=%.3f 4=%.3f\n\n", acc(m3a), acc(m3b), acc(m4)))## Accuracy 3A=0.576 3B=0.562 4=0.535cat("Model 3B (coef / odds ratio / p):\n")## Model 3B (coef / odds ratio / p):co <- summary(m3b)$coefficients
print(round(data.frame(Coef=co[,1], OddsRatio=exp(co[,1]), p=co[,4]), 4))## Coef OddsRatio p
## (Intercept) -0.1877 0.8289 0.7497
## Fscore 0.1733 1.1893 0.0000
## size -0.0349 0.9657 0.3545
## pb -0.0504 0.9509 0.1852
## roa 1.0593 2.8845 0.5176
## debt_r -0.2773 0.7578 0.29977 Multinomial Logistic Regression
library(nnet)
Dm <- D[!is.na(D$Cat5), ]; Dm$Cat5 <- factor(Dm$Cat5)
mm <- multinom(Cat5 ~ Fscore+roa+cfo_ta+d_roa+d_gm+d_ato, data=Dm, trace=FALSE)
cat("5-class accuracy =", round(mean(predict(mm)==Dm$Cat5),3), " (baseline 0.20)\n\n")## 5-class accuracy = 0.253 (baseline 0.20)cat("Relative Risk Ratios (vs class 1):\n"); print(round(exp(coef(mm)),3))## Relative Risk Ratios (vs class 1):## (Intercept) Fscore roa cfo_ta d_roa d_gm d_ato
## 2 0.617 1.097 0.262 0.058 0.391 1.521 0.838
## 3 0.279 1.274 0.037 0.006 1.313 0.986 0.097
## 4 0.184 1.388 0.334 0.007 11.197 0.924 0.111
## 5 0.267 1.297 14.706 0.018 0.615 1.006 0.409z <- coef(mm)/summary(mm)$standard.errors
cat("\np-values:\n"); print(round((1-pnorm(abs(z)))*2,4))##
## p-values:## (Intercept) Fscore roa cfo_ta d_roa d_gm d_ato
## 2 0.067 0.0797 0.6901 0.1533 0.6854 0.1525 0.9080
## 3 0.000 0.0000 0.3484 0.0147 0.9144 0.9627 0.1314
## 4 0.000 0.0000 0.7511 0.0162 0.3936 0.7977 0.1504
## 5 0.000 0.0000 0.3941 0.0388 0.8383 0.9856 0.54888 Visualizations
fit <- lm(fut_ret ~ Fscore, data=D); b<-coef(fit); r2<-summary(fit)$r.squared
plot(jitter(D$Fscore,1.5), D$fut_ret, pch=16, col=rgb(.17,.42,.72,.25),
xlab="Piotroski F-score", ylab="Next-quarter return (%)",
main="F-score vs. Future Stock Return", ylim=c(-70,80))
abline(fit, col="#c0392b", lwd=2.2)
legend("topleft", bty="n", text.col="#c0392b",
legend=sprintf("y = %.3f + %.3f x\nR2 = %.4f (p=%.3f)",
b[1],b[2],r2,summary(fit)$coefficients[2,4]))
f2<-lm(fut_ret~roa,data=D)
plot(D$roa, D$fut_ret, pch=16, col=rgb(.18,.52,.34,.25),
xlab="ROA (quarterly)", ylab="Next-quarter return (%)",
main="ROA vs. Future Stock Return", ylim=c(-70,80))
abline(f2, col="#c0392b", lwd=2.2)
legend("topleft", bty="n", text.col="#c0392b",
legend=sprintf("y = %.2f + %.2f x\nR2 = %.4f", coef(f2)[1],coef(f2)[2],summary(f2)$r.squared))
g <- tapply(D$fut_ret, D$Fscore, mean)
barplot(g, col="#3182ce", xlab="F-score", ylab="Mean next-quarter return (%)",
main="Average Future Return by F-score"); abline(h=0)
D$yqn <- D$year + (D$quarter-1)/4
tr <- tapply(D$Fscore, D$yqn, mean)
plot(as.numeric(names(tr)), tr, type="o", pch=16, col="#6b46c1",
xlab="Year", ylab="Mean F-score", main="Industry Average F-score Over Time")
write.csv(round(summary(m0)$coefficients,4), "out_OLS_simple.csv")
write.csv(round(summary(m1)$coefficients,4), "out_OLS_model1.csv")
write.csv(round(data.frame(OR=exp(co[,1]),p=co[,4]),4), "out_logit_3B.csv")
write.csv(round(exp(coef(mm)),4), "out_multinom_RRR.csv")
cat("Result tables exported.\n")## Result tables exported.